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Abstract 


Traditionally high-aspect ratio triangular/tetrahedral meshes are avoided by CFD re- 
searchers in the vicinity of a solid wall, as it is known to reduce the accuracy of gradient 
computations in those regions and also cause numerical instability. Although for certain 
complex geometries, the use of high-aspect ratio triangular/tetrahedral elements in the 
vicinity of a solid wall can be replaced by quadrilateral/prismatic elements, ability to 
use triangular/tetrahedral elements in such regions without any degradation in accuracy 
can be beneficial from a mesh generation point of view. The benefits also carry over to 
numerical frameworks such as the space-time conservation element and solution element 
(CESE), where triangular/tetrahedral elements are the mandatory building blocks. With 
the requirement of the CESE method in mind, a rigorous mathematical framework that 
clearly identifies the reason behind the difficulties in use of such high-aspect ratio triangu- 
lar /tetrahedral elements is presented here. As will be shown, it turns out that the degree 
of accuracy deterioration of gradient computation involving a triangular element is hinged 


on the value of its shape factor + = sin? a1 +sin? a2 +sin? a3, where ai, @2 and az are the 
internal angles of the element. In fact, it is shown that the degree of accuracy deterioration 
increases monotonically as the value of ~ decreases monotonically from its maximal value 
9/4 (attained by an equilateral triangle only) to a value < 1 (associated with a highly 
obtuse triangle). By taking advantage of the fact that a high-aspect ratio triangle is not 
necessarily highly obtuse, and in fact it can have a shape factor whose value is close to the 
maximal value 9/4, a potential solution to avoid accuracy deterioration of gradient compu- 
tation associated with a high-aspect ratio triangular grid is given. Also a brief discussion 
on the extension of the current mathematical framework to the tetrahedral-grid case along 
with some of the practical results of this extension is also provided. Furthermore, through 
the use of numerical simulations of practical viscous problems involving high-Reynolds 
number flows, the effectiveness of the gradient evaluation procedures within the CESE 
framework (that have their basis on the analysis presented here) to produce accurate and 
stable results on such high-aspect ratio meshes is also showcased. 


I. Introduction 


In the multi-dimensional space-time conservation element and solution element!~° (CESE) method, tri- 
angles and tetrahedral mesh elements turn out to be the most natural building blocks for 2D and 3D spatial 
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grids, respectively. As such, the CESE method is naturally compatible with the simplest 2D and 3D unstruc- 
tured grids and thus can be easily applied to solve problems with complex geometries. However, because 
(a) accurate solution of a high-Reynolds number flow field near a solid wall requires that the grid intervals 
along the direction normal to the wall be much finer than those in a direction parallel to the wall and, 
as such, the use of grid cells with extremely high aspect ratio (10° to 10°) may become mandatory, and 
(b) unlike quadrilateral/hexahedral grids, it is well-known that accuracy of gradient computations involving 
triangular/tetrahedral grids tends to deteriorate rapidly as cell aspect ratio increases, the use of triangular/ 
tetrahedral grid cells near a solid wall has long been deemed impractical by CFD researchers. ° 

In view of (a) the critical role played by triangular /tetrahedral grids in the CESE development, and (b) 
the importance of accurate resolution of high-Reynolds number flow field near a solid wall, a comprehensive 
and rigorous mathematical framework that clearly identifies the reasons behind the accuracy deterioration 
as described above has been developed for the 2D case involving triangular cells. By avoiding the pitfalls 
identified by the 2D framework, and its 3D extension, it has been shown numerically that the accuracy 
deterioration phenomenon identified above can indeed be overcome completely. 


II. A Summary of Key Theoretical Results 


Here, we provide a summary of the key practical results obtained from the current mathematical devel- 
opment. The mathematically-inclined reader is referred to Section III for the more rigorous mathematical 
treatment and details. 


Figure 1. General notations of the vertices, internal angles, lengths of the sides and altitudes of a triangle. 


To proceed, let (i) AP, P2P3, depicted in Fig.1, be a triangle lying on the x — y plane with its area, 
A(AP, P2P3) > 0 (2.1) 


and (ii) for each & = 1,2,3, let a, be the internal angle associated with the vertex, Py, I, be the length of 
the side facing the vertex P,, and hy be the length of the altitude originating from P;,. Then, (i) 


1 
(ii) 
I, > O and hy > 0,k = 1,2,3 (2.3) 
and (iii) 
T > Q1,Q2,a3 > 0 anda; +a9+03 =7 (2.4) 


Because of Eq.(2.3), the aspect ratio 7 of AP; P2P3 can be defined as 


}>0 (2.5) 
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i.e., 7 is the maximal value of the three positive ratios enclosed within the braces. In turn, by substituting 
a result of Eq.(2.2), ie., hy = 24, k&=1,2,3 into Eq.(2.5), one has 


ly 


ly Ip l 1 2 2 2 
(2A/l,)’ (2A/ly)’ (Ain) a xqmaxt (lr) ’ (lz) , (13) } >0 (2.6) 


7 = max{ 


Let (k1, k2,k3) be a permutation of (1, 2,3) such that In, > Ine > ly3. Then, with the aid of Eq.(2.3), we 
have 
lei = leo = les > 0 (2.7) 


As an example, in the case where lg > 1, > 13,k, = 2,k2 = 1, and k3 = 3. However in the case where 
1, = I, > Is, one has either (i) ky = 1,k2 = 2, and k3 = 3, or (ii) ky = 2,kg = 1 and kg = 3. With the aid of 
Eqs.(2.2) and (2.7), Eq.(2.6) now implies that 


 a* Maye. te 


= — = >0 2.8 
2A lev her Ae ee) 


Note that Eqs.(2.2) and (2.7) imply that lj; and hy, respectively, are the lengths of the largest 
side and the shortest altitude of AP, P2P3. According to Eq.(2.8), the ratio of these two lengths 
is the aspect ratio 7 of AP, P2P3. 

To recast 7 as a function of the internal angles a1,a2, and a3, note that, by observing Fig.1 and using 
Eqs.(2.3) and (2.4), one has the following relations: 


hy = lg sina3z = 13 sina > 0 
hg = Ig sina, = 1; sina3 > 0 (2.9) 


and hg = 1, sina2g = lg sina, > 0 
Next, by using Eq.(2.9) and the definition of k,,k2, and k3, one concludes that 

her = leo Sin agg = Ig sin agg > O (2.10) 
In turn, by substituting Eq.(2.10) into Eq.(2.8), one has 


lea bea 


” lp sin Ak3 Lx3 sin Ak2 ( ) 


Moreover, with the aid of a result of Eq.(2.4), ie., 
1> sina, > 0,k = 1,2,3 (2.12) 


Eq.(2.3) and the law of the sines imply that there exists a constant coefficient 3 > 0 such that 


lk l lz 
= kl =- k2 =— k3 = B > 0 (2.13) 
SIN Wk1 SIN Wrp2 SIN Wk3 
i.e., 
lp, = Bsinap, > 0, leo = Psinazeg > 0, and 1x3 = Bsinax3 > 0 (2.14) 
As a result of Eqs.(2.7) , (2.12) and (2.14), we have 
1> sinag > sinagg > sinaz3 > 0 (2.15) 


Moreover, by substituting Eq.(2.14) into Eq.(2.11), one can cast 7 as a function of a%1,%2 and az3, i.e., 


sin Wp1 


—— : (2.16) 
SIN Mp2 * SIN WKB 
As will be shown immediately, 7 can also be expressed as a function involving only ax2 and axz3. 
With the aid of Eq.(2.4) and the definition of k,,k2,and k3, we have 
T > Qi, AK2,%3 > O and ag, + Ogg + AK3 = 7 (2.17) 
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Thus, with the aid of some trigonometric identities, we have 
sin @y1 = sin (7 — Qp2 — Qx3) = Sin (@x2 + x3) = (Sin AxZ2)(COS ax3) + (Sin a3 )(COS Az2) (2.18) 
Substituting Eq.(2.18) into Eq.(2.16) and using Eq.(2.17), we have 
"= cotarge + cotargs3, O < axR2,QK3 < 7 and agg + aK3 < 7 (2.19) 


With the above preliminaries, the following proposition is proved in Appendix A. 
Proposition 1: Assuming Eq.(2.17), Eq.(2.15) implies that either 


T 
9 > Onl > an2 > ar3 > 0 (case Eq.(2.17)) (2.20) 


or 
T > Ap > and _ > 1 — Ap > An. > AR3 > 0 (case Eq.(2.17)) (2.21) 


must be true. In turn, Eqs.(2.16) and (2.19)—(2.21) imply that: (i) 


def 2 

n>" Fa (2.22) 

(ii) 
1 = min if and only if ay = ag = a3 = . (2.23) 

(iii) 
li = 2.24 
biyy 304 4 ae ( ) 

and (iv) 

n>>1 if and only if min{ay,a2,a3} =arn3 <1 (2.25) 


According to the above proposition, the aspect ratio 7 (i) attains its minimal value 2/./3 if and 
only if AP, P2P3 is equilateral, and (ii) attains a very large value if and only if the value of 
one of the internal angle of AP, P2P3 has a very small value. 

Next, let the “shape factor” of AP, P:P3 be defined by 


ay det in? a, + sin? ag + sin? a3 (2.26) 


It is shown in section III and Appendix B that (i) 


9 
0<7<5 (2.27) 
(ii) 
£9 
Y= Ymax Fi © a, = a2 = a3 3 26] =lg=l1, >0 (2.28) 
(iii) 
y 7 07 6 mar{ay,a2,03} 3 1 (2.29) 
and (iv) 
T 
max{ay, a2,a3} = 5 >7y=2 (2.30) 


Hereafter, (i) the symbol “@ ” is used to indicate that the statements given on its left and right sides are 
equivalent, while the symbol “ => ” is used to indicate that the statement given on its left side implies 
that given on the right side; (ii) as an example, y > 0* states that ~ approaches 0 from the range > 
0 while mar{a,,a2,a3} > m7 states that mar{a,,a2,a3} approaches 7 from the range < 7; (iii) as 
an example, the symbol {a1,a2,a3} denotes a set containing elements aj,a2, and a3, and, as a result, 
{a1, a2,a3} = {b1, b2, bs} indicates that the collection of elements a,,a2, and a3 is identical to that of the 
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elements 6), b2, and b3; and (iv) as an example, (a1, a2, a3) = {a1, {a1, a2}, {a1, a2, a3}} so that (a1, a2,a3) = 
(b1, ba, bs) © ay = by, k = 1,2,3. 

According to Eqs.(2.27)—(2.30), (i) y attains its maximal value 9/4 if and only if AP, P2P3 is equilateral; 
(ii) y approaches its minimal value 0 if and only if the largest internal angle of AP; P2P3 approaches 7~ (and 
therefore its two other internal angles approach 0*), e.g., the highly obtuse triangle case depicted in Fig.2(a); 
and (iii) 7 * 2 and thus has a value relatively close to maz if one of the internal angles of AP; PP3 is close 
to the right angle even if one of its internal angles approaches 0 such as the case depicted in Fig.2(b). Note 
that Ig/l3 = sinag/sina3 © ag/a3 if 0 < aez,a3 < 1. Thus, depending on the ratio ag/a3, the ratio of 
the lengths of the two sides P, P3 and P; P2 can take any value > 0 even for the highly obtuse triangle case 
depicted in Fig.2(a). 


P, 
pa 1 
P; Ps ee dt 
0 < a2, a3 «1 P; P; 
and 0<a, « 1, @, ~~ and 
0<n-a,<«1 a, =nm/2 


(a) (b) 


Figure 2. Representations of two kinds of high-aspect ratio cells: (a) highly obtuse, and (b) nearly right-angle triangle. 


Even though the two triangles depicted in Fig.2 have shape factors of vastly different values — one of them 
approaches the minimal value while another is relatively close to the maximal value of 9/4, both triangles 
have large aspect ratio. This is due to the fact that both meet the condition that min{a1,a2,a3} < 1 
which, according to Eq.(2.25) is the necessary and sufficient condition for 7 > 1. 

In Section ITI, the relation between accuracy deterioration of gradient computation and the shape factor 
y is established using a procedure briefly described below. First, for each k = 1,2,3, and at any given time 
level, let (i) 6, be the numerical value of a scalar function ¢ assigned to the vertex P;, of a given triangle 
AP, P2P3 lying on the x — y plane, and (ii) (az, yz) be the coordinates of the point P,. Second, let V4, 
the gradient vector of ¢ associated with AP, P:P3, be evaluated in terms of ¢,, 2%, yx, k = 1,2,3, using 
the standard linear interpolation procedure defined by Eqs. (3.2), (3.5), (3.6) and (3.21). Third, for each 
k = 1,2,3, let Ad, denote the numerical error of ¢, at the given time level introduced through a time- 
marching scheme while it is assumed that the coordinates (xx, yx), k = 1,2,3, be specified values with no 
numerical errors. Fourth, with the above assumptions it is shown in Eqs.(3.24), (3.35)—(3.37) and (3.40) 
that the numerical errors €1,€2 and e€3 of the directional derivatives along the side directions BP, ; PsP, ; 
and PP of AP; P2P3, respectively, are dependent only on [x and Adz, k = 1,2,3, while, according to Eqs. 
(3.26), (3.40)—(3.43), (3.46), (3.56), (3.57) and (3.66)—(3.69), the numerical errors Av, and Av, of 0¢/0x 
and 0¢/0y, respectively, are functions of Adz, x, and yx, k = 1,2,3 only. Fifth, the error norm for €1,€2 and 
€3 is defined as \/(e? + € + €%)/3 while that for Av, and Avy are defined as \/[(Av,)? + (Av,)?]/2. Sixth, 
as shown in Eq.(3.208), the error amplification factor R for gradient computation in turn is defined as the 
ratio between the error norm for Av, and Av, and that for €;,€2 and es. Seventh, according to Eq.(3.209) 
and other associated equations given in Section III, it turns out that the factor R is a function of the internal 
angles a1, @2 and as, and the errors €1,€2 and €3 only. In fact, according to Eqs.(3.167)—(3.170), the 2 x 2 
real symmetric matrix H, which appears in Eq.(3.209) is a function of a1,a2 and a3 only, while, according 
to Eqs.(3.54), (3.108), (3.112), (3.115) and (3.118), the 2 x 1 real column matrix w, which also appears in 
Eq.(3.209) is a function of a1,a2,a3,€2 and €3 (note: €1,€2,€3,Q1,Q@2 and ag are linked by Eq.(3.89), as 
such, given a set of a1,@2 and a3 , only two of €),€2 and €3 are independent parameters). Finally by using 
Eq.(3.209) and the Rayleigh-Ritz theorem(Ref.[8], p.43), one concludes that for, any combination of €1, €2 
and €3, (i) 


o-(y)SR< Vo+(y),0<7< 9/4 (2.31) 
where, 4 
mye 7 [3 + 2,/(9/4) — | 0<7< 9/4 (2.32) 
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are the two eigenvalues of H;, and (ii) each of the lower and upper bounds \/o_(7) and \/o+(7) in Eq.(2.31) 
can be attained by R by some special combinations of €1,€2 and €3. As such \/o_(y) and \/o+(7) are the 
greatest lower bound and the least upper bound of R, respectively. 

According to Eq.(2.32), both eigenvalues of H, are functions of the shape factor y only, even 
though Hj itself is a complicated matrix function of ai,a2 and ag defined by Eqs.(3.167)— 
(3.170). Also, by using Eqs.(2.27) and (2.32), it is shown in Section III that, as the value of 7 decreases 
from its maximal value 9/4 to its minimal limit 0+, (i) the value of o4(7) increases monoton- 
ically from 1 to +00; and (ii) the value of o_(-y) decreases monotonically from 1 toward the 
limit value of 1/2. Thus, with the aid of Eqs.(2.28) and (2.29), one has (i) 


o(y=lso(y=HlseMe=ah 4 the 2x 2 identity matrix 
9 a F (2.33) 
= 7 7 = QA] a2 a3 3 1 2 > 
(ii) 
1 ; 9 

5 <a_(y)<1l<o.(y), if0<7< 1 (2.34) 

(iii) 
1 
a(y) > Cy 6 04(7) 3 +00 6 ¥ 3 07 & mar{ay, 02,03} 3 1 (2.35) 


and (iv) 
3 
max{ay, a2,a3} = 7 S7y=25 Voy) = i: = 1.225 (2.36) 


With the aid of Eqs.(2.33)—(2.36), Eq.(2.31) implies that, (i) the error amplification R for gradient 
computation has the constant value 1 for any combination of €,,¢€2 and ¢3, if and only if AP, P:P3 is 
equilateral - as such the most accurate gradient computation occurs in the case where the triangular grid 
is built from equilateral triangles, (ii) for the case in which the grid is built from triangles similar to that 
depicted in Fig.2(b), the least upper bound of R is only slightly greater than 1 - as such accuracy deterioration 
of gradient computation is rather slight for this case, and (iii) because \/o(y) approaches +00 as y > 0*, 
accuracy deterioration of gradient computation progressively becomes worse when triangles used in grid 
construction become more obtuse and it would reach an unacceptable level if these triangles become highly 
obtuse, similar to the triangle depicted in Fig.2(a). Fortunately, the triangular grids suitable for two- 
dimensional viscous simulations near a solid-wall boundary can always be constructed using triangles similar 
to that depicted in Fig.2(b). As will be shown in Section IV, for grids containing triangles similar to that 
depicted in Fig.2(b), the gradient evaluation procedures within the CESE method can be effectively used to 
improve numerical stability while maintaining numerical accuracy. 

At this juncture, note that even for a grid built from nearly right-angle triangles similar to that depicted 
in Fig.2(b), and thus the least upper bound of R is close to 1, accuracy of gradient computation is still 
dependent on the magnitude of the error norm associated with €,,€2 and €3. As such, for simulation of a 
high-gradient boundary-layer flow near a solid-wall, accuracy of gradient computation can be enhanced if 
(i) the grid is aligned such that the shortest side P,P3 of AP, P2P3 depicted in Fig.2(b) is normal to the 
solid-wall, and (ii) the grid-interval size | P, P3| be reduced so that the magnitudes of the errors of the flow 


directional derivatives evaluated along the P, P3 direction be small enough for an accurate simulation. 

Moreover, the reader should be warned that accuracy deterioration of gradient computation associated 
with a grid built from highly obtuse triangles cannot be overcome though a simple “fix” described below. 
Given the highly obtuse AP, PP; depicted in Fig.2(a), let point P, be the perpendicular projection of P, 
on P,P3, i.e., Py lies on P23 P3 with P,P, being perpendicular to P;P3. Because 0 < ag,a3 « 1, Py must be 
an interior point on P2P3, 1.e., 


824 >0, $34 >0 and ly = 593 = 524+ 83.4 (2.37) 


where, 


def »35—35— def 
Ski ,ko = Sko,ki = | Pra Pies | =) (hy ~ Lk)? a (Yk _ Yko)* >0 


(2.38) 
ky Akg and ky, ko = 1,2,3,4 
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with the understanding that x, and yz, respectively, be x— and y— coordinates of point Py, k = 1, 2,3, 4. 
With the above definition of point P,, one can assign a value ¢4 of the scalar function ¢ to the point P, 
where 3 
gy S 44, 4+ 49, (2.39) 
52,3 §2,3 
ie., ga is the linear interpolation of ¢2 and ¢3 along P2P3. Because each of AP, P:P, and AP, P3P, is 
a right-angle triangle with its value of shape factor y close to the ideal maximal value 9/4, one may be 
tempted to conclude that accurate gradient computation over AP; P,P, (AP; P3P,) can be achieved using 
Le, Ye and dy, k = 1,2,4(k = 1,3,4). As will be explained below, because a large error introduced by the 
linear interpolation Eq.(2.39), this conclusion is misleading. 
Let P; denote a point in the x — y — ¢ space (see Fig.4) with the coordinates (24, y4,¢1). Then, with 


the aid of (a) Eq.(2.39), (b) the fact that P, is the perpendicular projection of P, on P2P3, (c) the fact 


that S(7/2) = Ip @ the2x2 identity matrix for the matrix function S(a) defined in Eq.(3.53), and (d) 


Eqs.(3.52) and (3.55), one can show that (a) P, is a point on the plane I defined by Eqs.(3.5) and (3.6), and 
(b) 


(Avy)? + (Avy)? = (e1)? + (€0)? (2.40) 


where (i) Av, and Av,, respectively, are the numerical errors of 0¢/0x and 0¢/Oy referred to earlier and 
evaluated using (2%, yr, Oe), k = 1, 2,3, (ii) «: is the numerical error of the directional derivative evaluated 


along P,P3 using xp, yx and dy,k = 2,3, and (iii) €g denotes the numerical error of the directional 


derivative evaluated along P,P, using zz, yz, and ¢%,k = 1,4. Then, according to the previous 
discussions, the error amplification factor R can become very large and thus the case 


(Av,)? + (Avy)? = (€1)? + (€0)? >> (1)? + (€2)? + (€3)? 


| (2.41) 
ie., (€9)” >> (€2)" + (€3) 


may occur for the highly-obtuse-triangle case depicted in Fig.2(a). As such the advantage of both AP, P2P, 
and AP, P3P, being right-angled triangles would be cancelled out by the fact that the error €) becomes 
unacceptably large. In fact, it can be shown that the simple fix described above has no impact on 
the accuracy of gradient computations other than round-off errors. 

Note that, as will be shown in a future paper, the mathematical development presented here for triangular 
grid can be extended in a straightforward manner to a tetrahedral-grid case, albeit the algebra becomes more 
complicated. In particular, the square of the error amplification factor R for a tetrahedral-grid case can still 
be cast into the form given on the extreme right side of Eq.(3.209) with the understanding that Hy is a 
3 x 3 real symmetric and positive-definite matrix [Ref.9, p.250] while w is a 3.x 1 real column matrix. 
In fact, it can be shown that, for a regular tetrahedron (i.e., a tetrahedron with all its four 
faces being equilateral triangles), H; is reduced to the 3 x 3 identity matrix, and as such 
the factor R = 1 for all possible combinations of the numerical errors associated with the 
directional derivatives evaluated along all six edge-directions of the tetrahedron. In other 
words, no accuracy deterioration of gradient computation occurs for a grid built from regular 
tetrahedrons. 

Moreover, it can be shown that, for a tetrahedron in which three right internal angles share 
a common vertex (trirectangular tetrahedron), the least upper bound of R is 2 ~ 1.414, a 
number slightly larger than 1. As such accuracy deterioration of gradient computation will be 
mild for a grid constructed from such tetrahedrons. 


Ill. A Mathematical Framework for Identifying the Cause and Cure of 
Accuracy Deterioration Associated with Triangular Grids with High 
Apsect Ratio 


As a preliminary consider Fig.3(a). Here points P,, P2, and P3 are the vertices of a triangle lying on the 
x—y plane with their coordinates, relative to the Cartesian x — y coordinate system shown in Fig.3(b), being 
(21, y1), (2, y2) , and (x3, y3) respectively. Let (i) 


A(AP,P2P3) = the area of the triangle AP; P,P; > 0 (3.1) 
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90° counterclockwise 
rotation 


P, 


(1) 


X 
(a) (b) 


Figure 3. General representation of (a) the triangle, along with its vertices, coordinates and internal angles; and (b) 
Cartesian coordinate system. 


(ii) the determinant of a square matrix M be denoted by |M|_ , (iii) 


ge te ed be 


3-21 Y3—- Yl 


= (x2 — £1) (ys — yi) — (@3 — £1) (ye — y1) (3.2) 


and (iv) the rotation from the positive direction of the x—axis to that of y—axis is in the counterclockwise 
direction with an angle of 90°.Then it can be shown that: (i) 


for the case in which the boundary of AP, P:P3 traced in the P; + P,; — P3 + P, direction forms a 
counterclockwise loop, and (ii) 
6 = —2- A(AP; P2P3) < 0 (3.4) 


for the case in which the boundary of AP, P)P3 traced in the P, + P, — P3; > P, direction forms a 
clockwise loop. 


P 


xX 
Figure 4. Representation of the x — y — ¢ space. 
Next, for each k = 1, 2,3, let (i) d, be a scalar value assigned to the spatial point P,, and (ii) P, denote 
a point in the x — y — ¢ space depicted in Fig.4 with the coordinates (x, yx, ¢,). In the following, first we 


will show that points P;, Pj, and P3 lie on the plane in the x — y — @ space represented by 


g=ax+byt+ec (3.5) 


where, 
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g2-or Ya-y t2-%, b2-1 TQ Yo 2 
3 — - 3-2 _ 
a lB 1 93 vay ae cca $3 ae ow [343% (3.6) 
6 6 6 
Proof: By definition, a plane in the x — y — ¢ space is represented by 
d\x + dgy+d3¢+d4,=0 (3.7) 
where dj, d2,d3, and d4 are real constants with 

(di, dz, ds) 4 (0, 0,0) (3.8) 


i.e., at least one of d;,d2, and d3 is non-zero. 
‘Next we will show that Piy Po, and Ps can only lie on a plane represented by Eq.(3.7) with ds 4 0. To 
prove by contradiction, let P,, P2, and P3 lie on a plane represented by Eq. (3.7) with ds = 0. Then 


dy, =F doy ial d4 => 0 (3.9) 

dix doy2 d4 =0 (3.10) 
and 

d\x3 doy3 d4 =0 (3.11) 


By subtracting Eq.(3.9) from Eqs.(3.10) and (3.11), respectively, one has 


(x2 — 21)d, + (y2 — yi)d2 = 0 


3.12 
and (a3 = x1)dy =e (ys == y1)d =0 ( ) 
On the other hand, according to Eqs.(3.2)—(3.4), we have 
fas ae eS (3.13) 
%—%1 Y3— YW 


In turn, according to elementary algebra, Eqs.(3.12) and (3.13) = d; = dy = 0. Because Eq.(3.7) does not 
represent a plane in the x — y — @ space if dy = dz = d3 = 0, points P,, Py, and P3 can only lie in a plane 
represented by Eq. (3.7) with d3 4 0. 


Let d3 #0. Then Eq.(3.7) = (i.e., is true if and only if) 
g=azrt+by+c (3.14) 


where q q d 
def 1 def 2 def 4 
a = ——, bl) = ——, andcd = —-— 


; 3.15 
ds 3 ds ( ) 


Note that, according to the above derivation of Eq.(3.14), the plane intersecting any three points P,, Pa, 
and P3 with A(AP, P2P3) # 0 (ie., Pi, Po, and P3 are not collinear points on the x — y plane) must be 
represented by Eq.(3.14) with the coefficients a’, b’, and c’ satisfying the conditions 


aa, +by +c = 91 
a'r + b'y2 +c = b2 (3.16) 
and a’ xy i b'y3 ar d = 3 


Because (i) Eq.(3.16) represents a system of three linear equations for the coefficients a’, b’, and c’, and (ii) 
with the aid of Eq.(3.13), we have 


rm yi il X Y1 I 
rq yo lj=|to-2, yo-y Of = 
x3 y3 1 3-2, y3—-y. O 


T—-7%1 Y2-Y1 
3-21 Y3— Y1 


=540 (3.17) 
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according to elementary theory of linear system of equations, the coefficients a’,b’, and c’ can be uniquely 
determined in terms of xz, yx, and dp, k = 1,2,3. In fact, with the aid of Eqs.(3.6) and (3.17), we have (i) 


di yi i P1 Yi 1 
g2 yo 1 g2—-b1 yo-y O d2—¢1 Yyo-y 
$3 y3 1 é3—-o1 ys-yi 0 d3—- 1 y3-Y 
= = = = 3.18 
= 5 5 5 = oe) 
(ii) 
xy gy 1 Ly or 1 
x2 g2 1 t2—-2, h2—-h, O 2-2, b2-h1 
rz $3 1 73-2, $3—¢, 0 3-21 63-1 
bo = = = =b a 
5 5 5 ot?) 
and (iii) 
Tm oy 1 
rT Yy2 2 
C= tee) =¢ (3.20) 
Thus the plane defined by Eqs.(3.5) and (3.6) is the same plane defined by Egqs.(3.14) and (3.18)—(3.20). 
QED. 


Hereafter let the plane defined by Eq.(3.5) and (3.6) be denoted by T. Then, on I’, Eq.(3.5) implies that 
the gradient vector 
= e O => O = i = 
Vé def sees + aa = aé; + bey (on T) (3.21) 
where é€;, and é, are the unit vectors in the c— and y— directions respectively. As such, on I’, (a) Voisa 
constant vector; and (b) both 0¢/0z and 0¢/0y are constant coefficients. 
Next, for each k = 1, 2,3, let = 
Ph = Grex + Yrey, K=1,2,3 (3.22) 


ie., P, is the position vector of point P, on the « — y plane. Also, for each (k,,k2) with ky # kg and 
ky, ke = 1,2,3 ; let 


> de > > 5 a 
Py, Pee = Pr, — Pry = (hp — ©, ee + (Yeo — Yer GG, hi # hy and ky, ko =1,2,3 (3.23) 
and 
e > 
Sk1,ko os Pr Pre = V (rk _ Tie 7 + (Yeo _ Yk, )?; ky # ko and ky, ko = 1, 2, 3 (3.24) 


— 
(Note here k; and kg are defined independent of Eq.(2.7)). According to Eq.(2.23), PPro is the displace- 
ment vector joining points P;,, and P,, and pointing towards Py, from Py, . Also by combining Eqs.(3.1) 
and (3.24), one concludes that 


Sky,ko = Sko,ky > 0, ky # kg and ky, kg = 1,2,3 (3.25) 
As a result, one can define 
a 
ef Pe, Pi e - 
Brisket and «ik, ke, = Pin Pha # ko and ky, ko = 1,2,3 (3.26) 
Ski ,ke $k1,k2 


Moreover, Eqs.(3.23)—(3.26) = (a) each €;, 4, is an unit vector, i-e., 


|Ek, ko | = 1, ky x ko and ky, ko = 1,2,3 (3.27) 


and (b) 
Eko ky = —Eky kg ONO [ky ,ky = —béky ko» Ki # ko and ki, ko = 1,2,3 (3.28) 
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As an example, consider the plane I which intersects points P,, P2, and P3. Let Vo be the constant 
gradient vector on T. Then Eqs.(3.21), and (3.23) > Vode Pr, Peo = Ok. — Ok, Which, with the aid of Eq. 
(3.26) in turn > 


Vb © Eby ky = Why kgs ki # ky and ky,kg =1,2,3 (on plane L) (3.29) 
In other words, for each pair of ky and kg with ky A ke and ky, kz = 1, 2,3, te, x, is the directional 
derivative of @ along the &,,,,, direction. The value of px,,~, can also be considered as the 


“slope” of P;,, Py, in the x — y — @ space. 
Moreover, Eqs.(3.25), (3.26) and (3.28) imply 


ee ae (3.30) 
$1,2 
— 13,2 = 12,3 = Pea (3.31) 
52,3 
and 
— 11,3 = 13,1 = os (3.32) 
53,1 
In turn, Eqs.(3.30)—(3.32) imply 
$1,2° M1,2 + $2,3° 2,3 + $3,1° 3,1 = 0 (3.33) 
and 
2,1 + M12 = 13,2 + H2,3 = 11,3 + 3,1 = 0 (3.34) 


Because the line segments P, Pz, P2P3 and P3P, are the three sides of AP; P;P3 , given the scalar values 
$1, G2 and os, /11,2, 2,1, 3,2; 2,3, #41,3 and 431 represent all possible directional derivatives of @ which can 
be evaluated along the spatial directions aligned with the three sides of AP, P2P3. Because the six directional 
derivatives are linked by four independent conditions given in Eqs.(3.33) and (3.34), any one of them can 
be determined in terms of any two independent directional derivatives among them. In fact, any two of 
them can be chosen as the independent directional derivatives, as long as they are evaluated 
along two different sides of AP, P2P3 and thus are not linked by one of the three conditions 
given in Eq.(3.34). 

For each & = 1, 2,3, let Ad; denote the numerical error of ¢, at some given time level introduced through 
a time-marching procedure. Moreover, (i) let the coordinates (x,,y,) of point Py,k = 1,2,3. be given 
fixed values that do not vary during the time-marching procedure, i.e., the numerical errors (Ax;,, Ay;,) of 
(tr, Ye), k = 1, 2,3, are assumed to be zero, and (ii) for any pair of ky and kz with ky 4 ko and ky, ko = 1, 2,3, 
let Apy,,~. Genote the error of 4x, ,~. introduced as a result of the errors Ady, and Ag, of o%, and dx, 
respectively. Then Eqs.(3.30)—(3.34) > 


Adz — A 
= Ape1 = Api.2 — ye Oe (3.35) 
$12 
Ads — A 
— Apsz,2 = Ape,3 = Ags — Ade (3.36) 
$2.3 
Ad, —A 
— Apis = Aus, = oa (3.37) 
53,1 
81,2 ° Apo + 82,3 - Ape3 + 53,1 -Apu31 = 0 (3.38) 
and 
Ape + Api2 = Aps,2 + Ape3 = Api3 + Asi =0 (3.39) 
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Because the six numerical errors Ap, .~.,k1 A ke and ky, ke = 1,2,3 are linked by four independent condi- 
tions given in Eqs.(3.38) and (3.39), any one of them can be determined in terms of any two independent 
errors among themselves. Any two of these errors can be chosen as independent errors as long as they 
are not linked by one of the three conditions given in Eq.(3.39), i.e., the numerical errors of any two 
directional derivatives evaluated along two different sides of AP, P2P3 can be chosen as the 
two independent numerical errors. 

To simplify the following developments, let 


€y at Ape3; €2 = Apsz,1 and €3 = Api,2 (3.40) 
and 


~ def ~ def , ~ def 
1 = €2,3, €2 = €3,1 and e€3 = €1,2 (3.41) 


Then, with the aid of Eqs.(3.28) and (3.39), Eqs.(3.40) and (3.41) > 


Aps3,2 = E15 A“i3 =—-€2 and Ape1 =—€3 (3.42) 


and 


€3.2 = —é), €13 = —€ and €21 = —€3 (3.43) 


At this juncture, note that the integers 1, 2, 3, appear in Eq.(3.40) as component-identification indices. 
Moreover under the index permutation, 


(1, 2,3) > (2,3,1) (3.44) 
(i.e., the integers 1, 2, and 3, respectively are replaced by 2, 3, and 1), Eq.(3.40) will turn into the form 


€9 ey Aps.1, €3 py Apti2 and €, = Apia.3 (3.40a) 


ie., with the exception of them being presented in different order, the definitions given in Eq.(3.40)a are 
identical to those given in Eq.(3.40). Hereafter, an equation or a set of equations that possesses such a 
property is said to be invariant under the index permutation of Eq.(3.44). In fact, each similar definition set 
to be given in Eqs.(3.109)—(3.111) is also invariant under the same index permutation. 

Note that, because the numerical errors (Av,, Ayx) of (ag, yx), = 1,2,3 are assumed to be zero in 
a time-marching procedure, according to Eqs.(3.22)—(3.26), for each pair of ky and kz with ky # kg and 
ky, kg = 1,2,3, the unit vector €;,,, is a constant vector with no numerical error during a time-marching 
procedure. As such with the aid of Eqs.(3.5), (3.21) and (3.40)—(3.43), and the definitions: (i) 


Vy, = “0 =a and vy, 2 se =b (onT) (3.45) 
and (ii) 
Ay © Av,é, +Av,é, (on) (3.46) 


where Av, and Av,, respectively, denote the numerical errors of vy, and vy. One concludes from Eq.(3.29) 
that 


Ape & = €x,k=1,2,3 (onT) (3.47) 
Also, because of the assumption Eq.(3.1) , (i) Eqs.(3.26) and (3.41) > 


Ek, x tk, if ky x ko and ky, ko = 1, 2, 3 (3.48) 


and (ii) the internal angles a1, v2 and a3 of AP, P2P3 depicted in Fig.3(a) must satisfy the condition (Note: 
the symbol “€” < “belongs to” ): 


(a1, a2, a3) €E Da (3.49) 


where 
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Da = {(a1, a2, a3)| 0 < a1,02,a3 <7 and a; + ag +a3 = 7} (3.50) 


Moreover, with the aid of Eq.(3.46), we have 


|x| SF Abe Av = y/(Avz)? + (Avy)? (3.51) 


and thus 


(|x|): = (Av,)? + (Avy)? (3.52) 


With the above preliminary, we will prove the following theorem: 
Theorem 1: Let (i) a1,@2, and ag be the internal angles of AP, P)P3 , shown in Fig.3(a) and thus 
(a1, A2, a3) € Da, (ii) 


ee 
S(a) # ( ; eal O<a<n (3.53) 
sin* a \ cosa 1 
and (iii) 
a2 (") (Qe (") and é; (") (3.54) 
€3 €1 €2 
Then 
2 
(|A>]) = (é,)'S(ax)é, for each k = 1,2,3 (3.55) 


where for each k = 1, 2,3, (é,) is the transpose of €,. 


4 


> 


Ck 


01. k =1, 2,3 
Xx 


Figure 5. Representation of a unit vector in x — y space. 


Proof: First note that, for each k = 1, 2,3, (i) & is the unit vector defined by Eqs.(3.26) and (3.41), and 
(ii) the angle 6, , shown in Fig.5, is uniquely defined by &, through the conditions: 


Ex = (cos 6% )Ez + (sin Oy )éy, k=1,2,3 (3.56) 
and 
T>O,>—m7, k=1,2,3 (3.57) 
Also it follows from Eq.(3.56) that 
Ek, © Eke = (cos Ox, )(cos O,,) + (sin 0%, )(sin Ox.) = cos(Ox%. — On), ki, ke = 1, 2,3 (3.58) 
Because €,,k = 1,2,3, are unit vectors, Eqs.(3.48) and (3.58) => 
(Ex, e Ek, ) = | cos(Oz., = Ox, )| <1 lif ky # ky and ky, ko = 1,2,3 (3.59) 


which, in turn, > 
sin(O,,. — O%,) #0 if ky A ky and ky, kg = 1,2,3 (3.60) 
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Note that, with the aid of Eq.(3.59), it can be shown that Eq.(3.48) <= Eq.(3.60). Moreover, because 


Eq.(3.57) > 
21 > Dk _ key > —2n, if ky x ko and ky, ke = 1,2,3 


As such, assuming Eq.(3.57), Eq.(3.48) <= Eq.(3.60) = 


Oby — On, A £m and Op, — On, #0 if ky # ky and ky, ko = 1,2,3 


Substituting Eqs.(3.46) and Eq.(3.56) into Eq.(3.47), one has 


(cos 6, )Av, + (sin@,)Av, = e, for any k = 1,2,3 


On the other hand, given any pair of k, and ka with ky 4 ko and ky, ky = 1, 2,3, in turn Eq.(3.63) > 


(cos 0%, Ave + (sin Op, Avy = €x, 


ky & ko and ky, ko = 1,2,3 
(cos O04, Ave + een, 1 # ke and ky, ke 52, 


Because, with the aid of Eq.(3.60) one has 


ROB Mis SUAVE tei nop \(doady, eda, ain Oa Seite, Oe 20 


cos 6, sin Oz, 


if ky x ko and ky, ke = 1,2,3 
Eq.(3.64) can be inverted to yield the following result 


6 = M(ky, ko)é(ki,ko) if ky # ky and ky, kg = 1,2,3 


where (i) 
7, AV, 
Avy 
(ii) 
def 1 sin 0; — sin 6; 
NCR er a 2 1) ky #ke and ky, ko =1,2,3 
(Fi, ka) sin(O;, — 9x, ) fo cos Ox, 1% ke ie 
and (iii) 


é(k1, ka) & (“) whi & ko and ky, ko = 1,2,3 
ko 


Note that: (i) Eqs.(3.54) and (3.69) > 
é, = &(2,3), é = é(3, 1) and é5 = (1,2) 


(ii) let [M(k1, k2)]’ be the transpose of M(ky, ke), then Eq.(3.68) > 


sin? (Oz, = Ox, ) = cos(Ox,, = Ox, ) 1 
for any (ki, kg) with ky x ko and ky, ko = 1, 2, 3 


They Ma RM eS ( 1 ~08() — i) 


and (iii) Eqs.(3.52), (3.66), (3.67) and (3.71) > 


(e|)” = (Ava)? + (Avy)? = (6)'5 = [ela he)]! M(B, a JM (hr, eC ha) 
= [é(k1, ko)|'T (ki, ka )é(k1, ka), ky x ko and ky, ko = 1; 2, 3 


where, 6* denotes the transpose of 6. 


NASA/TM—2017-219534 14 


(3.61) 


(3.62) 


(3.63) 


(3.64) 


(3.65) 


(3.66) 


(3.67) 


(3.68) 


(3.69) 


(3.70) 


(3.71) 


(3.72) 


Note that: (i) Eq.(3.69) > 


(ka, ky) = P é(ky, ka), for any (k1, ke) with ky # ko and ky, ko = 1,2,3 (3.73) 
where the 2 x 2 permutation matrix 
pee (3.74) 
1 0 
has the property 
PP =P (3.75) 


where P~! is the inverse of P; (ii) Eq.(3.71) > 
T(ko,k1) =T (ki, ka), br # by and ki, ke = 1,2,3 (3.76) 
(iii) with the aid of Eqs.(3.74) and (3.75), Eq.(3.71) > 
P*T (ky, ko)P =T (ki, ke), ha & kp and Ai, ka = 1,2,3 (3.77) 
and (iv) Eqs.(3.73), (3.76), and (3.77) > 


[é(ke, k1)|"T (ke, k1)é(ko, ki) = [€(k1, ka)]'T'(k1, ke) e(k1, ke) 


(3.78) 
ky # ko and ky, ke _ 1,2,3 


Because of the relation Eq. (3.78), Eq.(3.72) represents only three independent relations. As such, with the 
aid of Eq.(3.70), one concludes that Eq.(3.72) = 


(|AvP?) = (A:)'T2.3)& = (@)'TB, Ne = ()'TU, es (3.79) 


By comparing Eq.(3.79) with Eq.(3.55) it is seen that the proof of Theorem 1 is completed if one can show 
that 
T(2,3) = S(a1), T(3,1) = S(a2) and T(1, 2) = S(asz) (3.80) 


To prove Eq.(3.80), note that Eqs.(3.25), (3.26), and (3.41) > 


(3.81) 


In turn, with the aid of Fig.3(a) and Eq.(3.81) , one has 

ee een =) |= 

P, P3 e P,P |P.P5| [Pr P2| cos ax 
|PiPs| |PiP| |PLPS| PiP| 


On the other hand, by using Eqs.(3.58) and elementary trigonometry, one has 


—€2 @ €3 = 


= cos ay (3.82) 


— €2 @ €3 = — cos(63 — 62) = —cos(O2 — 43) (3.83) 
By using Eqs.(3.82) and (3.83) along with elementary trigonometry, one has 
— cos(62 — 63) = cosa, and sin?(62 — 63) = sin? ay (3.84) 


Moreover, by using Eqs.(3.60), (3.71) and (3.76), one concludes that (i) sin(02 — 63) 4 0, and (ii) 


sin? (O2 A 3) 1 


T(2,3) = T(3,2) = = ( - 1 — €08(82 — ) (3.85) 


Combining Eqs.(3.53), (3.84), and (3.85) one arrives at the first part of Eq.(3.80), ie., T(2,3) = S 
Moreover, because the new version of proof emerging from the above proof for T(2,3) = S(a1) 
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after it undergoes the index permutation of Eq.(3.44) remains valid, one will arrive at the second 
part of Eq.(3.80). Similarly one can prove the third part of Eq.(3.80). Because Eqs.(3.79) and (3.80) > 
Eq.(3.55), the proof of Theorem 1 is completed. QED. 

As a preliminary for the following developments, note that, by using (i) Fig.3(a), (ii) Eq.(3.25), and (iii) 
the law of the sines, one concludes that: (i) 


sina, > 0,k = 1,2,3, for each (a1, a2,a3) € Da (3.86) 


and (ii) there exists a parameter 6 > 0 such that 


$23 $31 — §1,2 
sin a1 sinag sinag3 


where Dj is defined in Eq.(3.50). Next, by using Eq.(3.40) and an equivalent of Eq.(3.87), i-e., 


= 6 > 0 for each (a1, a2,a3) € Da (3.87) 


$23 = Bsinay, $31 = Ssinaz and s12=sinag (6 > 0 and (a1, a2,a3) € Da) (3.88) 
Eq.(3.38) can be recast as 
(sin a1 )e, + (sin ag)eg + (sina3)e3 = 0, (a1, 02,03) € Da (3.89) 
Next, for a given (a1, Q2,a3) € Da, let 


A(ay1, @2,03) = {(€1, €2, €3)|€1, €2 and €3 be real parameters 


3.90 
satisfying Eq.(3.89)}, (a1,a2,a3) € Da oy) 


ie., for a given (@1,Q@2,Q@3) € Dag, A(a1,a2,a3) contains all the elements (€1,€2,€3) which 
satisfy Eq.(3.89). In turn, with the aid of Eq.(3.86), one concludes that, for any (€1, €2,€3) € A(a1, a2, a3) 
with (a@1,a2,a3) € Da, any one of €,,€2, and €3, can be uniquely determined in terms of the other two by 
using Eq.(3.89). In fact, for any (€1, €2, €3) defined from any given (€1, €2,€3) € A(a1, a2, a3) using Eq.(3.54), 
Eq.(3.89) < any one of the following six relations: 


€1 — J (a1, 2, 03 )é2, €9 = [J (a1, 02, 03)|~*e1 (a) 
€2 = J (a2, a3, 01) és, €3 = [J (a2, a3, A1)| 12 (b) (a1, Q2, 43) € Dy (3.91) 
and €3 = J (a3, 01, @2)é1, €] = [J(a3, a1, a2)|1é3 (c) 


Here, for any (a1, @2,a3) € Do [which } (a2, 03,01) € Da = (a3, 01, 2) € Dal, (i) 


ah ad os 
J (a1, W2, 03) = : ( ie >) ,(Q1, 02,03) € Da (3.92) 
sina2 \ sinag 0 
and (ii) [J(a1,a2,a3)|~+, is explicitly given by 
=i 1 0 sin a, 
[J(a1, a2, a3) = : . . ? (Q1, A2, a3) S Da (3.93) 
sina; \—sinag —sinag3 


As a preliminary for a key future development, given any (a1,Q2,a3) € Da, next we will prove the 
following relations: 


S(a1) 
S Q2) 
and S(a3) 


a1) = [(J(or, 02,03) J’ S(a2)(J(a1,a2,a3))-? (a) 
4)" $(a3)(J(a2,03,01))7? (b) > (a1, 02,03) € Da (3.94) 
4 S(ax)(J(a3,01,02))7! (0) 


Proof for Eq.(3.94): By using Eqs.(3.53) and (3.93), one has 


(a1) = ))- 
(a2) = [(J(a@2,03,01))~ 
= ye 


[(J(as, 1, 2 


1 


7) ey) 
sin a1 -sin* ag 


[(J(a1, a2, a3))~4]* S(a2)(J(a1, a2, 03))7! = 
(3.95) 


sin? a sin a2(sin a3 — sin a, cos a2) 
sin a2(sina3 — sina, cosa2) sin? a; + sina3(sina3 — 2sin a, cos az) 


’ (Q1, @2, a3) € Da 
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To simplify Eq.(3.95) note that Eqs.(3.49) and (3.50) = a3 = 7 — aj — a2 and thus 
sin a3 = sin(m — ay — a2) = sin(a, + a2) = sina, cos a2 + sin a2 cos a (3.96) 
In turn, by using Eq.(3.96) and some trignometric identities, one has 
sin a2(sin a3 — sin ay Cos a2) = (sin? a2) cos ay (3.97) 


and 
sin? a; + sina3(sina3 — 2sina; cosaz) = sin? ag (3.98) 


Next, by substituting Eqs.(3.97) and (3.98) into Eq.(3.95) and then using Eq.(3.53) again, one arrives at 
Eq.(3.94)(a). Because Eq.(3.94)(b) is the image of Eq.(3.94)(a) under the index permutation of Eq.(3.44), by 
carrying out this permutation over the proof for Eq.(3.94)(a), one has the proof for Eq.(3.94)(b). Similarly 
one can prove Eq.(3.94)(c). QED. 

Note that an immediate result of Eq.(3.55) is 


(€)'S(a1)e = (€2)'S(a2)é = (é3)'S(a3)é3, (a1, 2,03) € Da (3.99) 


In the following, Eq.(3.99) is proved directly using Eqs.(3.91) and (3.94). 
Proof for Eq.(3.99): By using the second part of Eq.(3.91)(a), one has 


(€2)'S(ag)éo _= (é,) [(J(ar, O92, a3))~*] : S(a2)(J(ar, Q2, a3))*é, (3.100) 
In turn, by substituting Eq.(3.94)(a) into Eq.(3.100), one has 
(€2)'S(az)ée => (€1)'S(ay)éy (3.101) 


ie., the validity of the first equality sign in Eq.(3.99) has been proved. Similary, one can prove the validity 
of the second equality sign. QED. 

Moreover, Eq.(3.91) implies that any (€1,€2,€3) which are defined in terms of the same (€1,€2,€3) € 
A(a1, @2,03) using Eq.(3.54) has the property that €, = 0 for any k = 1,2,3 @ €; = €2 = €3 = O. As 
such, any (€1, €2,€3) defined using the same (€1, €2,€3) € A(a1,@2,a3) must meet one of the following two 
mutually exclusive conditions: 


(a) €, =0, for each k = 1,2,3 (3.102) 


and 
(b) é, #0, ie., (€,)°€, > 0, for each k = 1,2,3 (3.103) 


With the aid of Eq.(3.54), it is seen that Eq.(3.102) = 
€& =—S€&gHeg = 0 (3.104) 


According to Eqs.(3.35)—(3.37) and (3.40), Eq.(3.104) implies that the numerical errors of the directional 
derivatives of ¢ evaluated along the three sides of AP, P2P3 are all zero. 
Because (i) Eq.(3.102) <= Eq.(3.104), (ii) Eq.(3.104) and the condition 


e& © ea)? + (e2)? + (es)? > 0 (3.105) 


eh 


are mutually exclusive, and (iii) any (€1, €2, €3) defined in terms of the same (€1, €2,€3) € A(a1, a2, a3) must 
belong to one of the two exclusive cases, Eqs.(3.102) and (3.103), one concludes that Eq.(3.103) <= Eq.(3.105). 
Because the case Eq.(3.104) is numerically unrealistic and, according to Eqs.(3.54) and (3.55), it leads 
to the trivial result Av, = Av, = 0 (which, according to Eqs.(3.40)—(3.43) and (3.45)—(3.47), simply states 
that the numerical error of the vector Vo on the plane [ is zero if the numerical errors €1, €2 and €3 of the 
directional derivatives of ¢ evaluated along the three sides of AP; P2 P3 are all zero), hereafter, for any given 
(a1, Q2,a3) € Da, unless specified otherwise, we will consider only the case Eq.(3.103), i.e., only the case 


(€1, €2, €3) € A'(a1, @2,@3), (a1, a2, a3) € De (3.106) 
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where, 
A’ (a1, 02, 3) ae {(€1, €2, €3)|(€1, €2,€3) € A(a1, a2, a3) and also satisfies 


(3.107) 
the condition Eq.(3.103) or its equivalent «* > 0} for any (a1,a2,a3) € Da 
Next, with the aid of Eqs.(3.86), (3.89) and (3.105), we will prove the following theorem: 
Theorem 2: Let (a) (a1,a2,a3) € Da, and thus Eq.(3.86) > 
O,, ff =“ > 0, k Al and kl =1,2,3 (3.108) 
sin a, 
(b) 
—_— 2 . 
B, 2 (1 Ga Ca Cis (3.109) 
Crea-Cig 14+ (C13) 
—_— 2 . 
pele See (3.110) 
Ca3°Co1 1+ (Ca) 
and 
2 . 
ae + (C31)"  C3,1 C32 (3.111) 
2 31°C32 1+ (C32) 
(c) 
dip 21+ (Cia)? + (Cia)? - 21 (3.112) 
dot 2 1+ (Ca3)? + (Cor)? Ao 21 (3.113) 
at def def 
Ag4 = 1+ (C31)? + (C32), As- = 1 (3.114) 
(d) 
w, = : VAi4 Cia VA1+ C13 (3.115) 
Vv (C1,2)? + (C13) C13 —C1,2 
W, 2 : VAz+ C23 VAz+ C21 (3.116) 
V (C2,3)? + (C2,1) Co. —C2,3 
and 
Ww, VA3+ C31 A384 C3,2 (3.117) 
V(C3,1)? + (C3,2) C3,2 —C31 
and (e) 
de = Weer, b= 1,2,3 (3.118) 


where €,,& = 1,2,3, are defined by Eq.(3.54) and by assumption, belong to the case Eq.(3.103). Then we 
have: 
(a) The real symmetric matrices E,, £2 and E3 are all positive definite [Ref.9, p.250]. 


(b) 


(e*)? = (€,)' Exe, > 0 for each k = 1,2,3 and each (€1, €2,€3) € A’(ay, a2, a3) (3.119) 
(c) 
Ex = (We)*We, k = 1,2,3 (3.120) 
and (d) 
(e*)* = (W)'vx > 0 for each k = 1,2,3, and each (€1, €2,€3) € A’(a1, 2, 03) (3.121) 


Proof: Note that, for each k = 1,2,3,A,4 and A,— defined in one of Eqs.(3.112)—(3.114) are the eigen- 
values of the real symmetric matrix E;,, defined in one of Eqs.(3.109)—(3.111). On the other hand, according 
to Eqs.(3.108) and (3.112)—(3.114), we have 


Art > Ar—- = 1, for each k = 1, 2,3 (3.122) 
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Thus, for each k = 1, 2,3, the eigenvalues of FE; are positive. As a result, by definition [Ref.9, p.250], for each 
k = 1,2,3, the real symmetric matrix FE, is positive definite. In turn, this coupled with Eq.(3.103) implies 
that [Ref.9, p.250]: 


(€,)' Exe, > 0 for each k = 1,2,3, and each (€1, €2,€3) € A’(a4, a2, 03) (3.123) 
Next, note that Eqs.(3.89) and (3.108) > 
€, = —(Ci2 €2 + C13 €3), with C12 >0 and C)3 > 0 (3.124) 
Then, by using Eqs.(3.54), (3.105), (3.109), (3.123) and (3.124), one concludes that 
(e*)? = [1 + (C12)*](€2)? + [1 + (C1,3)7](€s)? + 2C1,2 » C1, - (€2es) 


ee (3.125) 
= (€,)"E1€, > 0 for each (€1, €2, €3) € A’(a1, a2, a3) 


ie., the & = 1 case of Eq.(3.119) has been proved. To prove the k = 1 case of Eqs.(3.120) and (3.121), note 
that (i) 


e 1 . de 1 
a,2 ee ee C13 (3.126) 
Vv (C1,2)? + (Ci,3)? \C1,3 J (C12)? + (C1,3)? \-Ci2 
are the eigenvectors of E, with eigenvalues 4,4 and ,_ respectively; and (ii) 
(24) = (G_)ar_ =1 (3.127) 


Moreover, as expected from Eq.(3.122), and a matrix theorem, i.e., two eigenvectors of a real symmetric 
matrix with different eigenvalues must be orthogonal to each other [Ref.9, p.222], we have 


(€14)'€1_ = (41_)"G14 =0 (3.128) 


Because €;4 and €,_ satisfy Eqs. (3.127) and (3.128), by definition [Ref.9, p.122], they form a pair of 
orthonormal eigenvectors of the matrix F. 


Let, 
US : ae 9 (3.129) 
V(C1,2)? + (Ci,3)? (C13 —C1,2 


i.e., the eigenvector €;, and €,_ of FE, respectively, are the first and second column of matrix U;. By using 
orthorormal property of €;4 and €,_, it can be shown that 


(U,)*U, = U,(U;)* = (U,)? = In & the 2 x 2 identity matrix (3.130) 


Thus Uj is a real symmetric orthogonal matrix [Ref.9, p.126] satisfying the relation 
(i) +=) =" (3.131) 


Moreover, because U; is formed using two orthonormal eigenvectors €;4 and €;— of FE), the latter can be 
diagonalized through a similarity transformation involving U, [Ref.9, p.223], i-e., 


7 a 0 Au «0 
U,) (EU, = = 3.132 
( 1) 141 ( 0 ) ( 0 ) ( ) 


where A; and A;_ are the eigenvalues of E, defined in Eq.(3.112). Moreover, because 414 > Ai- = 1 and 
thus \/A14 > 1, Eqs.(3.131) and (3.132) > 


E, =U, & ) (U,)-' = (U,)* & Mt : U; (3.133) 


1 


Moreover, because Eqs.(3.115) and (3.129) > 


Wi = (° — ) Uy (3.134) 
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Eq.(3.133) now implies the k = 1 case of Eq.(3.120), ie., 
E, =(u1) (\ = ) ( : U, = (Wi)'W, (3.135) 


Next, by combining Eqs. (3.125), (3.135) and (3.118), one has 
(e*)? = (€1)’Exé: = (€1)*(Wi)*Wiés = (Wiér)'Wieér = (1)*h1 > 0 for the case Eq.(3.103) (3.136) 


i.e, the k=1 case of Eq.(3.121) has been proved. 

To prove Eqs.(3.119)—(3.121) for the & = 2 case, note that new version of proof emerging from the above 
proof for the k=1 case after its component-identifications indices 1, 2, and 3 undergo the index permutation 
of Eq.(3.44) remain valid. In particular, after this permutation, Eqs.(3.125), (3.135) and (3.136), respectively, 
become oe 

(e*)? = (€9)’ Egéy > 0, Ey = (W2)'(Wo) and (e*)? = (2)*H2 > 0 (S137) 


ie., the k=2, case of Eqs.(3.119)—(3.121) has been proved. Similarly, the k=3 case can be proved by carrying 
out the index permutation of Eq.(3.44) over the proof for k=2 case. As such, Theorem 2 has been proved. 
QED. 

Next, the following corollary follows directly from Eqs.(3.86), (3.108) and (3.112)—(3.118). 
Corollary to Theorem 2: Let (a ,a2,a3) € Da. Then, (a) 


jokes iene Si eS (3.138) 
sin” a, sin ap 
where 
(a1, A2, a3) def sin? a, + sin” ag + sin? a3 > 0, (a4, 02,03) € Dy (3.139) 


is referred to as the shape factor of AP; P2P3, hereafter; (b) 


|W;,| “ the determinant of W, = —/ny = — VY ec 1k=1,2,3 (3.140) 


sin A, 


Thus, for each k = 1, 2,3, (W,)~! exists; (c) Eq.(3.118) = 


é, = (Wi) ~ de, k = 1,2,3 (3.141) 
and therefore 
Eq.(3.103) © oy F (") , Le., (We), > 0, for each k = 1,2,3 (3.142) 
and (d) 
1 ae C13 
(W,)' = a a ee . (3.143) 
VJ (C1,2)? + (C13) rama 
C2,3 
C. 
1 2,1 
(W2)? = ; 2| (3.144) 
V (C2,3)? + (C2,1) en —C2,3 
and 
1 ne C3, 
(W3)7! = ; oo (3.145) 
J (C3,1)? + (C3,2) cers —C31 
As a preliminary for a key future development, note that, by using Eqs.(3.93), 
(3.108), (3.115)-(3.117), (3.138)-(3.140) and (3.143)—(3.145), one has 
Gy 2 Ws (J (a2, 03,01))~"(Wa) 
= 1 ( sin a2 - sina3 —,/ysinay (3.146) 
(sin? a2 + sin? a1) (sin? a3 + sin? a1) V7sin ay — sin ag - sin a3 
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Go W; (Flos, 01,02))-(Wa)? 


= 1 (- sin a3: sina, —,/7sin a2 (3.147) 
(sin? a3 + sin” ag)(sin? a, + sin? a2) V7sin ag — sin ag - sin ay 
and 
def = = 
G3 = We (J(a1, @2,a3)) (Wi) * 
i ( sina, -sinag  —,/ysinag3 (3.148) 
(ein? a1 + sin? a3) (sin? ag + sin? a3) V7sin. a3 — sin q - sin a2 
Hereafter, only cases with (a1, a2,a3) € Da will be considered. 
Next, with the aid of Eqs.(3.146)—(3.148) and the identities: 
(sin? a2 + sin? a1) (sin? a3 + sin? a1) = sin? a2 - sin? a3 + ysin? ay (3.149) 
(sin? ag + sin? a2)(sin? a1 + sin? a2) = sin? a3 - sin? a, + ysin? a2 (3.150) 
and 
(sin? a; + sin? a3)(sin? ag + sin? a3) = sin? ay - sin* ag + ysin? a3 (3.151) 
One has 
(Ge)'(Ge) = Ia, k= 1,2,3 (3.152) 


Thus, for each k = 1,2,3,G, is a real orthogonal matrix [Ref.9, p.126], ie., 
(Ge) + =(Ge)’, &=1,2,3 (3.153) 
Next, by using Eqs.(3.91)(a)—(3.91)(c), (3.118), (3.141) and (3.146)—(3.148) one can show that 


de = Gsthi, 3 = Gitte and ti = Goy)3 (3.154) 


Proof for Eq.(3.154): By using Eqs.(3.118), (3.91)(a), (3.141) and (3.148) one has 


ho = Woe = Wo(J (a1, 02,03))~1é, = Wo(J (a1, a2,03))~1(Wi)7 dh = Gah 


ie., the first part of Eq.(3.154) has been proved. Similary, by using Eqs.(3.118), (3.91)(b), (3.141), (3.146) 
and (3.147), one can prove the second and third parts of Eq.(3.154). QED. 
Note that, by using Eq.(3.153), it can be shown that Eq.(3.154) = 


th = (Gs) "G2 = (Gs)"vho, do = (Gi) ~'b3 = (Gi) !b3 and Ws = (G2)7 hr = (G2)'hi (3.155) 
Also, by using Eqs.(3.106), (3.107) and (3.142), one concludes that, the current basic assumption that 
(€1,€2,€3) € A’(a1, a2,03) & 


dp # (") , ie, (We)'de > 0 for each k = 1,2,3 and each (a1, 02,03) € Da (3.156) 


At this juncture, also note that a result of Eq.(3.121) is 
(v1)'dr = (Go)! = (s)"ds (3.157) 


which can also be proved directly using Eqs.(3.153) and (3.154). Recall that Eqs.(3.153) and (3.154) are 
derived using (i) the relations given in Eqs.(3.91)(a)—(3.91)(c) which are derived from Eq.(3.89), (ii) the 
definition of Eq.(3.118) and its inverse relation Eq.(3.141), and (iii) the definitions G1,G2 and G3 given in 
Egs.(3.146)—(3.148). In other words, Eq.(3.157) represents a set of identity relations which follow directly 
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from Eq.(3.89). As such, with the aid of the fact that (€1, €2,€3) € A’(a1,a2,a3) & Eq.(3.156), Eq.(3.121) 
can be replaced by a simpler equivalent form, i.e., 


(2) Oe eth 2 (' 


) , Le., (d1)'d > 0 (3.158) 


Next, by using Eqs.(3.55) and (3.141), one has 
7 
~s 


5 (\Wel)? (Wi) da]? S (an) (We) de = (eb) Hah, k = 1,2,3 (3.159) 
where 3 
H,& 5((We) I'S (ax)(We) 1, & = 1,2,3 and (a1, 02,03) € Da (3.160) 


According to Eq.(3.53), for each k = 1,2,3,S(a,) is a real symmetric matrix. In turn, Eq.(3.160) = for 
each k = 1,2,3, Hz is also a real symmetric matrix. As such, for each k = 1, 2,3, the eigenvalues of H;, are 
all real [Ref.9, p.222]. Moreover, as will be shown, the matrices H,,H2 and H3 are similar [Ref.9, p.232], 
i.e., they are related by the following similarity transformations: 


Ay = (G3)-! Ay G3, H2 = (Gi)-! A3 G, and H3 = (Go)! Ay, Go (3.161) 


where G),G2 and G3 are the matrices defined in Eqs.(3.146)—(3.148), respectively. According to a matrix 
theorem [Ref.9, p.232], similar matrices such as H,, Hj and H3 have the same eigenvalues with the same 
multiplicities. In the following, first we prove Eq.(3.161). 

Proof of Eq.(3.161): By using Eqs.(3.160), (3.94)(a), (3.148) and (3.153) one has 


Ay = 5[((Wi)*]'S(a1)(Wi)* = S[(Wi)*]'[(J (a1, a2, 03))*]*S (a2)(F(a1, a2, 03))~"(Wi)* 


Do] w dol w 


(Ma) (F(a, 2, 4))-"]! Wa) [Way (@2)(Wa)- Wal T(an,a2,09))""Wa* 3 gy 


= [Wo F(a1, 02, 03))-"(Wi)*]¢5 [(W2)1]'S (2) (We) [Wa F(a, 02, 08))-()~) 
= (G3)'H2G3 = (G3) HaG3 


iLe., the first part of Eq.(3.161) has been proved. Next, with the aid of Eqs.(3.160), (3.94)(b), (3.146) and 
(3.153), it can be shown that the new version of Eq.(3.162) after it undergoes the index permutation of 
Eq.(3.44) is also valid, ie., the second part of Eq.(3.161) is also valid. Similarly, once can also prove the last 
part. QED. 

Obviously, Eq.(3.161) can be cast in the following equivalent form: 


Hz = G3H,(G3)~!, H3 = Gy H2(G1)~! and H, = G2H3(G2)~* (3.163) 
Moreover, by using Eqs.(3.153), (3.154) and (3.163), one can show that 
(di)! Hid = (v2)! Howe = (Ws)! Has (3.164) 
In turn, with the aid of Eq.(3.164), Eq.(3.159) can be replaced by a simpler form 
3,=> “ n 
5 (Vl) = (v1) ity (3.165) 


To study the eigenvalues of Hj, note that, with the aid of Eqs.(3.108) and (3.138), Eq.(3.143) > 


: sin a2 sin a3 
(W,)-! = daca vy — sinaa \ (a1 ,02,03) € D (3.166) 
1 a) a) sin a3 sinaa |’ , 13 oi : 
sin” Q@2 + sin” ag “VT sinay 


In turn, by substituting Eq.(3.53), (3.139) and (3.166) into Eq.(3.160), one has 


Mao 1 ©) fe aceeD, (3.167) 
2\e fi 
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where, 
def Sin? a2 + sin? ag + 2(cos a1) (sin a2)(sin a3) 


dy = 


Q1,A2,a3) € D 
(sin? a2 + sin? a3)(sin? ay + sin? a2 + sin? a3)’ (1,02, 03) ° 


def COS Q1 (sin? a3 — sin” Q2) 
ey = 


» (a1,42,03) € Da 


sin ay (sin? a2 + sin? as) (sin? a + sin? a2 + sin? a3) 


and 
def Sin” a2 + sin? ag — 2(cos a1) (sin a2) (sin a3) 


fi 


» (Q1,A2,a3) € D 

sin? a; (sin? a2 + sin? a3) (1,2, 43) . 

By definition, an eigenvalue o of H, is a root of the equation: 
the determinant of the matrix (H; — oI) =0 


Substituting Eq.(3.167) into Eq.(3.171), one has 
23 9 2 
oO — 5 (di + fide + Glas - fi — (ex) ]=0 


Thus, 
O = 044 OF 0 = O4_~ 


where 


oi4 = aa; + fit V(di + fi)? — Aldi - fr — (e1)?]} 


Note that, because 
(di + fi)? — 4[di - fa — (e1)7] = (di — fr)? +. 4(e1)? > 0 


(3.168) 


(3.169) 


(3.170) 


(3.171) 


(3.172) 


(3.173) 


(3.174) 


(3.175) 


the expression under the radical sign in Eq.(3.174) is non-negative, Thus, for any (a@1,a@2,a3) € Da, the 
eigne values 0,4 and o,_ of Hy are real always, a fact consistent with the fact that the eigenvalues of the 


real symmetric matrix H, must be real always [Ref.9, p. 222]. 
Next, with the aid of Eqs.(3.139) and (3.168)—(3.170), one has (i) 


dy: fi - (e1)? = 


(sin? ag + sin” a3)? — 4(cos? a1)(sin? a2)(sin? a3) — (cos? a1) (sin? a3 — sin? a2) 


(sin? a;)(sin? ag + sin? a3)?(sin? a; + sin® ag + sin” a3) 


(sin? a2 + sin? a3)? — (cos? a;)(sin? a2 + sin? a3)? 


(sin? a1) (sin? ag + sin? a3)?(sin? ay + sin? a2 + sin? a3) 


(sin? a1) (sin? ag + sin? a3)? 


(sin? a1) (sin? ag + sin? a3)?(sin? ay + sin? a2 + sin? a3) 
1 1 
= = Q1,a2,a3) € D 
(sin? a; + sin? ag +sin2a3)  Y’ (01, 02, a3) = 


and (ii) 
1 
7(sin? a1) (sin? a2 + sin? a3) 


d+fi= 


where. 
qi & (sin? a1) [sin? ag + sin? a3 + 2(cos a1)(sin a2)(sin a3)| 
+ (sin? a1 + sin? a2 + sin? a3)[sin? ag + sin? a3 — 2(cos a1)(sin a2)(sin a3) 
= (sin? a1) (sin? ag + sin? a3) + 2(sin? a1)(cos a1) (sin a2)(sin a3)+ 
(sin? a1 + sin? ag + sin? a3)(sin? ag + sin? a3) — 2(sin? a1) (cos a1)(sin a2) (sin a3) 
— 2(sin? ag + sin? a3)(cos a1) (sin a2)(sin a3) 
= (sin? a2 + sin” a3) [2 sin? a; + sin? a2 + sin? a3 — 2(cos a1)(sin a2)(sin a3) 
To simplify Eq.(3.178), next we will prove the identity: 


sin? ag + sin? a3 — 2(cos a1)(sin ag) (sina3) = sin? a1, (a1,a2,03) € Da 
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(3.176) 


(3.177) 


(3.178) 


(3.179) 


Proof of Eq.(3.179): By using Eq.(3.24) and Fig.3(a), the law of cosines can be expressed as 


(81,3)? + (81,2)? — 2(cos a1) 81,3 - 81,2 = (82,3), (a1, @2,03) € Do (3.180) 


Because $3.1 = $1,3 and o > 0, Eq.(3.179) now follows directly from Eqs.(3.88) and (3.180). Eq.(3.179) can 
also be proved directly without invoking the law of cosines. Because, (a1, a2,a3) € Da => a, = 7 —Q2— 43, 
we have 


sin? a, = sin? (a2 + a3) = (sin az cos a3 + sin a3 cos a2)” 


= (sin? a2) (cos? a3) + (sin? a3) (cos? az) + 2(sin a2)(sin a3)(cos a2) (cos a3) 
= (sin? ag)(1 — sin? a3) + (sin? a3)(1 — sin? a2) + 2(sin a2)(sin a3)(cos a2) (cos a3) 


= sin? a2 + sin? a3 + 2(sin a2)(sin a3) cos(ag + a3) = sin? ag + sin? a3 — 2(cos a1) (sin a2)(sin a3) 
i.e., Eq.(3.179) has been proved. QED. With the aid of Eq.(3.179), Eq.(3.178) now implies that 
q1 = 3(sin? a1) (sin? a2 + sin? a3) (3.181) 


In turn, by substituting Eq.(3.181) into Eq.(3.177), one has 


3 
dit+fi=—, (a1,@2,03) € Da (3.182) 


2 


Next, with the aid of Eqs.(3.176) and (3.182), Eqs.(3.174) and (3.175) imply that, for any given 
(a1, a2, a3) € Da, (i) 


ee ere ra £5 1G nl, Gavenon en. (3.183) 
and (ii) 
SG) = (= fy)? + Alex)? 20, (01,02, 08) € Ds (3.184) 
In turn, Eqs.(3.139) and (3.184) > 
0<7< : (a1, Q2,a3) € Da (3.185) 


Note that, it was shown earlier that H,, H2 and Hz are similar and thus they must have the 
same eigenvalues. As such, Eq.(3.183) implies that o1(y) and o_(y) must be invariant under 
the index permutation of Eq.(3.44) — A fact that leads to the expectation that they must be 
a function of parameters (such as ) which are invariant under this permutation. 

Moreover, it is shown in Appendix B that, by using Eq.(3.139), one can prove directly that, (i) Eq.(3.185) 
represents the range of y over Dg; (ii) 


9 
Y= 7 1 = 2 = a3 if (a1, 2,03) € Da (3.186) 


i.e., the shape factor 7 reaches its maximal value $ if and only if AP; P)P3 is equilateral; and (iii) for the 
case (@1,Q2,a3) € Da, 


> 07 & min{ay,a2,a3} > OF 
‘ aia (3.187) 
and max{a,, a2,03} > 1 


ie., for the case (a1, 02,03) € Da, y > 0F & the value of the largest internal angle approaches 7~ while 
those of the other two approach 0*. 

Moreover, by using Eqs.(3.183), (3.186) and (3.167)—(3.170), one can show that, for the case (a1, a2, a3) € 
Day G) 


Hi, = I for the case (a1, a2,a3) € Da (3.188) 
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(ii) 


o4(5) =land lim oi(y) =+00 (3.189) 


0+ 


(iii) 


9 
is undefined at y = Z 


do+(y) 
a _ ee] 9 (3.190) 


if0<y<- 
8 OA) VO/4) —7 1 
and (iv) 
; 9 
is undefined at y = — 
do_(7) ‘ 
ates (3.191) 
dy) _3[9=271-6VO/D=7) sey 2 
8 (97) V (9/4) — ¥ 4 
Because, 
9 9 
9 — 27 + 69/4) — y= 45 — 7) + 274+ 6/ (9/4) -y>0,0<7< i (3.192) 
Eq.(3.190) > 
do+(7) 9 
iy <0 if0<y¥<7 (3.193) 
On the other hand, because 
j-d26 (0 =9= (9 — 27 — 6/ (9/4) — 7)(9 — 27 + 6 (9/4) — 7) 
9 — 27 + 6,/ (9/4) — y (3.194) 
ti en) eee 
(9 — 27 + 6/(9/4) — 7) 4 
Eq.(3.191) => 
do-(y) _,; 9 
pu poito<y<} (3.195) 


To study the behavior of o_(y) in the limit of y + 0*, note that (i) 


li —2 A)-—y)= li — 1 
as (9/4) —y) =0 and Be 0 (3.196) 
and (ii) 
d|3 — 2 4) — 
B-2/9/)-7_ 1 ge 4. Gee (3.197) 
dy (9/4) — 7 oe . 


Then, with the aid of Eqs.(3.196) and (3.197), an application of L’Hospital rule over the function o_(¥) 
defined in Eq.(3.183) => 


3 1 3.2 #1 
li (y= > li = = 3.198 
ae +g, ( a 1*3 73 oan 


By using Eqs.(3.189)—(3.191), (3.193), (3.195) and (3.198), one concludes that: as the value of 7 decreases 
from ? to 0* , (i) the value of o4(7) increases monotonically from 1 to +00; and (ii) the value of o_ (7) 
decreases monotonically from 1 towards the limit value 1/2. Thus, one has (i) 


o(yy=lso(y=leyv=- (3.199) 
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(ii) 
1 9 
3 < 9-1) <1 <4) f0<y¥<q (3.200) 
and (iii) 
1 
a_(y) > Gy &a4(y) 7 +00 & y 3 OF (3.201) 
Note that Eq.(3.139) > 


y = 2 if mar{aj, a2,a3} = . and (a1, a2,a3) € Da (3.202) 
Thus, with the aid of Eq.(3.183), one concludes that 


3 


o+(y) = 04(2) = ; and o_(y) = o_(2) = ri 


if (3.203) 
if max{ay, a2,a3} = . and (a1,@2,a3) € Da 


At this juncture, note that Eqs.(3.139), (3.183), (3.185)—(3.187) and (3.199)—(3.201) imply that: (i) 
TW OT 

_ EN ae (eV ae -~\)=1] 204 

)=o(5, 5,5) = 08) =04(8) (3.204) 

(ii) 


1 
3 < 71- (01, a2, 03) = o-(7) <1 < o14(a1, a2, a3) = +(7) 


ee . (3.205) 
if (a1, @2, a3) al (a a? a) (ie., is Y< =) and (a1, a2, 3) € Da 
3°3°3 A 
and (iii) 
o1(a1,.02,08) = o-(9) > (5)” 014 (a1, 02,08) = 04 (9) 4 $00. 
(3.206) 


7 2 07 & min{ayz, a2,a3} + OF and maz{ay1,a2,a3} > 7, 

if (a1, Q2, 03) Ee Da 

However, by using the results presented above and the fact that an n x n real symmetric matrix possesses 

a set of n linearly independent real eigenvectors associated with its exclusive real eigenvalues [Ref.8, p.306], 


one concludes that, for each (a1, 2,03) € Da, there exist two linearly independent real eigenvectors _ (7) 
and y,4(7) such that 


Aydhi—-(y) =0-(y)i-(y) and Aida (y) = (yd () 


9 (3.207) 
for each y withO<y< Z and each (a1, @2,a3) € Da 


With the above preliminaries, we are ready to tackle the central question that motivates the current 
study. First, note that, by using Eqs.(3.51) and (3.105)-(3.107) 


act [3|Av| | (Ave)? + (Ary)?I/2 
_ Va o~ ne + (ea)? (es)3/B = Or (e126) € A'(, 02, 088) (3.208) 


By definition, R is the square root of the ratio of the two simple averages, i.e., the simple average of (Av,)? 
and (Av,)?, and that of (€1)*, (€2)? and (e3)?. As such R is a measure of the relative magnitudes 
of the error norms e*/\/3 and |Av|/V2. According to Eqs.(3.35)—(3.37), (3.40) and (3.105), e*/VW3 
is an error norm associated with the numerical errors of the directional derivatives of the 
scalar function @ evaluated along the three sides of AP, P2P3. On the other hand, according to 
Eqs.(3.45),(3.46) and (3.51), |Av|/V2 is an error norm associated with the numerical error of the 
constant gradient vector ¢ on the plane [ in the « — y— @ space, defined by Eqs.(3.5) and (3.6). 
As such R is a measure of how large the error norm associated with the gradient vector @ on 
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T is amplified from that associated with the directional derivatives of @ evaluated along the 
three sides of AP, P2P3. Thus, by its definition, R represents an error amplification factor. 

Numerical experiments reveal that the value of R could (but not necessary) become very large if AP; P2P3 
has a large aspect ratio. In the following, we will provide a mathematically rigorous explanation of this 
pitfall and also a way to avoid it even in a case in which use of a triangular grid with a large aspect ratio is 
unavoidable. 

To proceed, note that it was shown earlier that the assumption (€1, €2,€3) € A’(a1, a2, a3) = Eq.(3.156). 
Thus, by combining Eqs.(3.158), (3.165) and (3.208), one has 


— A a 
3/2)(\Av|)? _ (Wi) Aida OV. paves 
= ( = for any dy + ( ie.,(t1)*1 > 0) (3.209) 
(e*)? (hi )*hr 0 
In turn, with the aid of (i) Eq.(3.209), (ii) the fact that H, is a real symmetric matrix, and (iii) Rayleigh-Ritz 
Theorem [Ref.8, p.431], one concludes that, given any (a1, a@2,a3) € Da 


o_(7) < R? <04(7) for any y # (") (3.210) 
Because, Eqs.(3.199), (3.200) and (3.208) > 
i “ 0 
3< a_(y) <1 <o0,(4) and R>0 for any y 4 (") and any (a1,Q2,a3) € Da (3.211) 
Eq.(3.210) © 
o(¥) <R< Vox(y) for any 4 F (") and any (a1, 02,a3) € Da (3.212) 


Moreover, Eq.(3.211) and the Rayleigh-Ritz theorem also = for any (a1, @2,a3) € Da, 


(i) R= v/o_(7) for any yy with Hyd, = o_(y)hy and y, 4 (") [(a1, a2, a3) € Da] (3.213) 


and 


(ii) R= /o4(y) for any dy with Aid = 04(y)hr and yy # (") [(a1, @2,a3) € Da] (3.214) 


To study the implication of Eqs.(3.212)—(3.214), first consider the special case in which a1 = ag = a3 = §, 
ie., AP, P2P3 is equilateral. According to Eqs.(3.188) and (3.199), for this case, one has 
9 9 9 
y= we-(q) - o+(7) =landM=[2 (a, =a, =a03 = *) (3.215) 
As such, for this special case, (i) Eq.(3.212) > 


R=1 for any py 4 (") (@1 = ag =a3 = 5) (3.216) 
and (ii) 
Hy = Ini = hy for any hy A (") (a1 = a2 =a3 = =) (3.217) 


ie., any yy x (") is an eigenvector of H, with unit eigenvalue if AP; P,P; is equilateral. 
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Next consider the case in which 


(a1, Q2, a3) # ( ) and (a1, a2,a3) € Da (3.218) 


T 7 7 
3°37 3 
According to Eqs.(3.185) and (3.188), Eq.(3.218) = 


9 
0<7y< ri (3.219) 


Because (i) Eq.(3.200) implies that o_(y) < o+(7) for each y with 0 < 7 < 7, (ii) two eigenvectors of 
Hy associated with two distinct eigenvalues, respectively, must be linearly independent [Ref.8, p.268], and 
(iii) the 2 x 2 matrix H, can have at most two linearly independent eigenvectors, one concludes that, for 
each 7 with 0 < 7 < 3, the eigenspace of Hj associated with either of the eigenvalues o_(y) and o+(7¥) is 
one-dimensional. As such, for a given 7 with 0 < 7 < 3, (i) 1, which is 2 x 1 nonzero real column matrix 
by definition, is an eigenvector of H; associated with the eigenvalue o_(y) if and only if 


: : 9 
vy =r_wy_(7), r- £0 for any y with0 <7 < i (3.220) 


where r_ is any real number # 0, and %,_ (7) is the real eigenvector of H, introduced in Eq.(3.207), and 


(ii) W1 is the real eigenvector of H) associated with the eigenvalue o,(c) if and only if 


: : 9 
Y1 =re vir(y), r+ £0 for any y with 0<y7< Fi (3.221) 


where r, is any real number ¥ 0, and ~#14(7) is the real eigenvector of H, introduced in Eq.(3.207). In turn, 
by combining Eqs.(3.212)—(3.214) with the observations given in the above items (i) and (ii), one concludes 
that, for any y with0O << 3 (ii) the lower bound ,/o_(7) is attained by R if and only if #, is in 
the form specified by Eqs.(3.220) and (ii) the upper bound of \/o+(y7) is attained by R if and 
only if 1 is in the form specified by Eqs.(3.221). As such, for the case Eq.(3.218) (or equivalently the 
case Eq.(3.219)), Eq.(3.212) implies that \/o,(7) and \/a_(7), respectively, are the least upper bound and 
the greatest lower bound of R for each y with 0 < 7 < 2. 

As noted earlier, by its definition Eq.(3.208), R is a measure of how large the numerical error norm 
associated with the gradient vector ¢ on the plane [ is amplified from that associated with the directional 
derivatives of @ evaluated along the three sides of AP; P2P3. On the other hand, for each given y with 
0<y< $ Eq.(3.212) states that the value of R can vary between the greatest lower bound \/o_ (7) and the 


least upper bound \/o+(7). To minimize the possible value of the numerical error norm |Av|/V2 
associated with the gradient vector of ¢ on I, against a given error norm e* /+/3 associated with 
the directional derivatives of ¢ evaluated along the three sides of AP, P2P3, it is imperative 
to minimize the least upper bound \/a+(7). 

According to the discussion given preceeding Eqs.(3.199)—(3.201), within its domain 0 < y < 2, the value 
of o+(7) increases monotonically from 1 to +00 as 7 decreases from 7 to 0*. As such \/o+(7) (i) reaches 
its minimal value 1 if and only of y = 2 and (ii) approaches +00 if and only if in the limit of y — OT. 

Moreover, according to Eq.(3.186), y = 9 if and only if AP; P2P3 is equilateral. Thus ,/o+(7) reaches 
its minimal value 1 if and only if AP; P2P3 is equilateral, a result one would expect intutively. On the other 
hand, according to Eq.(3.187), y — 0* and therefore \/o,(y) approaches +00, if and only if the values of 
the two smaller internal angles of AP, P2P3 are both much less than one, as the case depicted in Fig.2(a). 
Luckily, even though AP; P;P3 depicted there must have a high aspect ratio, by no means the values of 
two of its internal angles must both be much less than one for a triangle with high aspect ratio. In fact, 
according to a discussion given in Sec.2, AP; P2P3 has a high aspect ratio if and only if the value of one of 
its internal angle is much less than one. As such, a triangle with the value of only one of its internal angle 
being much less than one, such as the case depicted in Fig.2(b), is also associated with a high aspect ratio, 
In fact, according to Eq.(3.202), y = 2 for any right triangle even if the value of its smallest angle is much 
less than one and thus it has a high aspect ratio. Moreover, according to Eq.(3.203), for this case, 


the least bound of R = \/o+(2) = iE = 1.22( for a right triangle). (3.222) 
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In other words, the least upper bound of R can be quiet close to its minimal value 1 even if the right triangle 
AP, P2P3 has a high aspect ratio. Thus, as long as one of its internal angles is close to the right angle, the 
least upper bound of the error amplification factor R associated with AP, P:P3 is quiet close to the minimal 
value 1 even if the triangle has a large aspect ratio. 

From the above theoretical discussions, one concludes that the accuracy of gradient calculation of a 
solution variable over a triangular spatial grid is tied to the shape factor y of individual triangles from which 
the grid is built. Specifically, if the values of y associated with the component triangles are close to its 
maximal value 2 (which is exactly attained if and only if the triangle is equilateral), the upper bounds of the 
error amplification factors R associated with the component triangles will be relatively close to its minimal 
value 1. It is also shown that the value of y associated with a triangle with a large aspect ratio can be 
close to the ideal maximal value 2 as long as one of its internal angles is close to the right triangle (ie., a 
90 degree angle). As such, accuracy deterioration of gradient calculation due to the high aspect 
ratio associated with the triangular grid used can be avoided if the grid is built from triangles 
with each of them possessing the property that one of its internal angles is close to the right 
angle. As will be shown in the following sections, the above theoretical predictions have been confirmed 
numerically. 


IV. Results and Discussion 


Numerical algorithms to circumvent numerical instability and accuracy issues for high aspect ratio tri- 
angular/tetrahedral meshes are being implemented in the software framework. The in-house CFD code 
framework ez4/d° based on various schemes that belong to the CESE family has been continuously under 
development since 2004. This software framework has been written using a combined object-oriented and 
generic programming paradigms in the C++ programming language. Light-weight object-oriented hierarchy 
is used in conjunction with heavy use of template classes and functions to allow compile time polymorphism. 
Different conservation laws can be plugged in with templates that represent physics. Currently, the software 
supports either triangular/tetrahedral or quadrilateral/hexahedral unstructured meshes. Both multi-thread 
(based on low-level POSIX thread) and message passing interface (MPI) paradigms are used to facilitate 
large-scale parallel computations. Each MPI process within a computational node can be executed in multi- 
thread mode to further enhance parallel performance, especially for a memory bound multi-domain layout. 
A communication map is used for data transfer among interface zones. For a large mesh, of the order of a 
billion elements, each unstructured block can be built with its own connectivity and nodes. A global commu- 
nication map is then used to join all the independent blocks in the parallel computations. This arrangement 
allows the grid generation process to always have a low memory requirement. The interfaces among blocks 
can be continuous or discontinuous. A continuous interface mesh ensures better solution accuracy for un- 
steady flow computations. Both second- and fourth-order CESE numerical schemes are implemented for 
general conservation laws including Euler and Navier-Stokes equations in the software framework. The time 
accurate local time-stepping (TALTS) scheme,’ ~° that allows all the elements in a mesh to march at an 
approximately uniform CFL number without violating flux conservation in time, is used to enhance parallel 
performance (smaller size elements do more computations than the larger ones) and numerical accuracy. 

High aspect ratio meshes are used routinely in viscous flow computations. In this section, several practical 
test examples are discussed to demonstrate the effectiveness of the CESE method to deal with high-aspect 
ratio meshes inside the boundary layer without any degradation in accuracy or numerical stability. The 
gradient evaluation procedures that enables CESE to use high-aspect ratio meshes are detailed in a previous 
paper.’ Furthermore, all the results featured here were computed with the use of the TALTS scheme, 
implemented in ez4d, to maintain computational efficiency and accuracy at the same time. 


IV.A. Acoustic Wave Propagation along High-Aspect Ratio Boundary-Layer Mesh 


The TALTS method is designed to handle large grid size disparities in a discretized domain. One frequently 
encountered problem in flow-physics is to simulate acoustic waves propagating through a viscous mesh. 
Computing acoustic waves using such a non-uniform mesh is quite challenging. Figure (6) shows a rectangular 
domain with a clustered triangular mesh with an aspect ratio (7) around 225 near the bottom boundary. A 
constant time step computation results in significant phase error due to the large CFL number disparities. 
The wave amplitude is also incorrectly damped. In contrast, the TALTS computation maintains a rather 
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Figure 6. Acoustic waves propagating through a domain with clustered mesh, showing density contours at a time 
instant. 


uniform CFL number (around 0.8) and thus results in much improved accuracy in both amplitude and phase 
throughout the domain. The results indicate that high aspect ratio triangular mesh poses no numerical issues 
for acoustic wave propagation. Capability to resolve acoustic or other physical waves inside the boundary 
layer is essential for direct numerical simulations of transitional or turbulent boundary layers. 


IV.B. Mach 3 Isothermal Laminar Boundary Layer 


Surface heating prediction has been one of the most important topics in many supersonic or hypersonic 
computational investigations. Accurately predicting the surface heat flux in a highly cooled boundary layer 
has important implications in vehicle thermal shield and aerodynamic design. Of interest is the isothermal 
boundary layer development behind the leading-edge bow shock of a blunt body. To isolate the boundary 
layer from shocks for numerical accuracy studies, a Mach 3 flow (with a free-stream temperature of 200 
K) over a highly cooled boundary layer with Tyau/Twatt,adiabatic = 0.2 is numerically computed. Due to 
the presence of the boundary layer, there is a weak leading edge shock that would make the free-stream 
conditions of the boundary layer deviate slightly from the Mach 3 conditions. The computational domain 
is a 1.2m x 0.1m rectangle (the height is about 5 boundary-layer thicknesses at the exit). A structured 
quadrilateral mesh with different sizes is sliced to form triangular elements for computations. The computed 
velocity and temperature profiles (with a mesh of 128,000 triangular elements and a maximum aspect ratio 
(n) of about 1500) at the location of Re = 10° are compared with compressible similarity boundary layer 
solutions in Fig.7, along with the surface heat flux distribution along the streamwise direction. The results 
show good agreement for laminar surface heat flux predictions. 


IV.C. RANS Computations for a Mach 2 Adiabatic Boundary Layer 


Meshes with high aspect ratio elements are most frequently used for Reynolds Averaged Navier-Stokes 
(RANS) computations. To resolve the viscous sublayer, log layer, and all the way up to the boundary layer 
edge, the near wall y* value has to be kept around or smaller than one. This step size is several orders 
of magnitude smaller than that along the strewamwise direction, which implies a very large aspect ratio. 
RANS computations are carried out using both Sparlart-Allmaras and Mentor’s SST models.'° Both skin 
friction coefficients and turbulent mean streamwise profiles are compared with results from the widely used 
NASA CFL3D code. A triangular mesh with about 105,000 elements and an aspect ratio (7) of around 3000 
obtained by slicing the quad mesh downloaded from the turbulent resource website,!° was used for a Mach 
2 adiabatic boundary layer. The agreement with structured mesh solutions is quite good (see Fig.8). It 
indicates the numerical dissipation treatment for the high aspect ratio unstructured meshes is adequate in 
resolving RANS equations. 


IV.D. Mach 11 Double Cone Benchmark Problem 


The CUBRC wind tunnel tests!! for hypersonic flows over a 15° — 25° double cone is often used to vali- 
date surface heating and the capability of CFD codes in predicting laminar or turbulent shock boundary 
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Figure 7. Comparison of Navier-Stokes solutions with compressible boundary-layer solutions at Re = 10° for a Mach 


3 flate-plate flow with Tya11/Twatl,adiabatic = 0.2 : (a) streamwise velocity (non-dimensionalized by u..) (b) temperature 
(non-dimensionalized by T,) and (c) heat flux distribution along the streamwise direction. 
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Figure 8. Comparison of RANS solutions with results from CFL3D for a Mach 2 flat-plate turbulent boundary layer 
flow with adiabatic walls: (a) skin friction coefficient (b) turbulent velocity profiles at Reg = 10,000. 
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layer interaction. Laminar calculations for the Run 35 case in Ref.[11] are carried out and the computed 
Mach number contours are shown in Fig.9 below. In all of the plots shown, the lengths have been non- 
dimensionalized with the length L of the first section of the cone. A triangular mesh with about 313,000 
elements and an aspect ratio (7) of about 1400 near the wall was used for the computations. The separation 
bubble around the corner and the big subsonic bubble around the second ramp are clearly captured. The 
cone surface pressure coefficient and Stanton number are compared with the experimental data in Fig.10. 
The agreement of Cp is quite good, while the current results under-predict the surface heat flux. A grid 
sensitivity study is currently underway to sort out the heat flux discrepancies. 
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Figure 10. Comparison of computed surface pressure coffeicent and Stanton number with experimental data for the 
Mach 11.3 flow over a double cone: (a) C, (b) Stanton number. 


V. Conclusion 


This work discusses the fundamental accuracy issues when a mesh with very large aspect ratio elements is 
used in CFD simulations, especially inside the viscous boundary layer. Sources of inaccuracy as an outcome 
of triangular shapes are identified theoretically, followed by a discussion on potential remedies. The current 
status about the application of the CESE method in steady or unsteady computations for acoustic waves 
or flow over a viscous boundary layer with large aspect ratio triangular meshes is briefly discussed, to 
demonstrate the effectiveness of the gradient evaluation procedures to deal with mesh-related geometrical 


difficulties. 
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Appendix A 


Proposition 1 presented in Section II will be proved here. To prove Eqs.(2.20) and (2.21), first note that 
Eq.(2.17) implies that each (ax, @%2, 3) Must meet one of the following two mutually exclusive conditions: 


(a) > max{ap1, M2, k3} > 0 with {ky, ke, kg} = {1, 2,3} (A.1) 


and 
(b) mt > max{agi, AK2, AK3} > 5 with {k1, ko, k3} = {1, 2,3} (A.2) 


Because, in the interval 0 < a < 7/2,sina@ increases monotonically from 0 to 1 as a increases 
from 0 to 7/2, it is seen that, for the case Eq.(A.1), Eq.(2.20) is a direct result of Eqs.(A.1) and (2.15). 
In the following, for the case Eq.(A.2), we will prove the Eq.(2.21) is also a result of Eqs.(2.15) and (2.17). 

First we will prove that, for the case Eq.(A.2), 


Apa A Max{Ag1, AK2, Ae3} and Ag3 A Mar{Ag1, Ak2,AK3} ( case Eq.(A.2)) (A.3) 
To prove the first part of Eq.(A.3) by contradiction, assume 
T > Apg = Max{ap1, K2, K3} > and case Eq.(A.2) (A.4) 


and let ne 
Opp = 1 — OK (A.5) 


Because Eq.(A.5) > apg = 7 — Ge2, Eq.(A.4) and the second part of Eq.(2.17), respectively, imply that 


> az > 0 for case Eq.(A.4) (A.6) 
and 
Ak2 = Aki + WKB (A.7) 
In turn, Eqs.(A.6) and (A.7) > 
5 > ap + a3 > 0 for case Eq.(A.4) (A.8) 


On the other hand, according to the first part of Eq.(2.17), agi > 0 and az3 > 0. Thus, Eq.(A.8) implies 
that 
> aki > 0 and ; > apg > 0 for case Eq.(A.4) (A.9) 


T 
2 
Next, because sin(z — a) = sina for any real number a, Eq.(A.5) implies that 
sin Qp2 = sin Wa (A.10) 
In turn, with the aid of Eq.(A.10), Eq.(2.15) implies that 
1> sinag > sin agg > sinaz3 > 0 (A.11) 


On the other hand, because, in the interval 0 < a < 7/2, sina increases monotonically from 0 to 1 as a 
increases from 0 to 7/2, for the current case Eqs.(A.6) and (A.9), Eq.(A.11) © 


. > an >On > ans >0 (case Eqs.(A.6) and (A.9)) (A.12) 
By substituting Eq.(A.7) into a result of Eq.(A.12), we have 
Qn. > OK2 = Ani +axz3 >0 (case Eqs.(A.6) and (A.9)) (A.13) 


which leads to the result 0 > ax3 and thus a contradiction to the first part of Eq.(2.17). As such the 
assumption Eq.(A.4) is false and therefore the proposition given in the first part of Eq.(A.3) is 
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required for consistency to Eqs.(2.15) and (2.17). Similarly, it can be shown that the proposition given 
in the second part of Eq.(A.3) is also required for consistency to Eqs.(2.15) and (2.17). 
Because of Eq.(A.3), for the case Eq.(A.2), consistency to Eqs.(2.15) and (2.17) requires that 


T > Ap = Max{ag, AK2, UK3} > 7 for case Eq.(A.2) (A.14) 


Because 
7 7 ____ def 
T> Oh 2 5 5 2 OR = 7 — Oe > O (A.15) 


one can show that, for the case Eq.(A.14), Eq.(2.17) © 


TT 

3 > Ak = T — Ak = Ang + AZ > Age > O 

and for case Eq.(A.14) (A.16) 
TT 

a > Ok = 7 — Ag = ARQ + KZ > AZ > O 


To search for the extra conditions on a@%1,@%2 and az3 required for consistency to Eq.(2.15), note that, 
by using the definition of Gz given above, and the relation sin(a — a) = sina for any real number a, it is 
seen that Eq.(2.15) = 

1> sindg > sinazg > sinagg > 0 (& Eq.(2.15)) (A.17) 


On the other hand, because (i) in the interval 0 < a < 4,sina increases monotonically from 0 to 1 as a 
increases from 0 to $, and (ii) Eq.(A.16) > 5 > agi > 0,5 > ax2 > O and § > ax3 > 0, for the case 


2 
Eq.(A.16), Eq.(A.17) & 
> Oe = 7 — Og > AK. > aKZ > 0 case Eq.(A.16) (A.18) 


As such, for the case Eq.(A.14), Eqs.(A.16) and (A.18) = Eqs.(2.15) and (2.17). On the other hand, it can 
be shown that Eqs.(A.16) and (A.18) = 


> — Op = On. + OK > Ane > AZ >0 case Eq.(A.14) (A.19) 


As such, for the case Eq.(A.14), it has been shown that Eqs.(2.15) and (2.17) = Eq.(A.19). Because, it 
was shown earlier that, the case Eq.(A.14) is the only subcase of Eq.(A.2) which can be consistent with 
Eqs.(2.15) and (2.17), one concludes that, for the case Eq.(A.2), Eqs.(2.15) and (2.17) = Eq.(A.19). 
Because, Eqs.(A.19) < (i) Eq.(2.21) and (ii) the second part of Eq.(2.17), it has been proved that, for the 
case Eq.(A.2), Eq.(2.21) is a result of Eqs.(2.15) and (2.17). 

It was shown earlier that, (i) as a result of Eq.(2.17), each (Qx1, x2, 43) Must meet one of the two 
mutually exclusive conditions Eqs.(A.1) and (A.2); and (ii) for the case Eq.(A.1), Eq.(2.20) is a result of 
Eq.(2.15). By combining the results with the result presented in the last paragraph, one concludes that 
Eq.(2.15) and (2.17) imply either Eq.(2.20) or Eq.(2.21). In the following, Eqs.(2.22)—(2.25) will be 
proved using Eqs.(2.16) and (2.19)—(2.21). 

Note that: (i) Eqs.(2.17) and (2.20) = 


5 > Qp1 > Ogg > 0 and 5 > az; >0 case Eqs.(2.17) and (2.20) (A.20) 
and (ii) Eqs.(2.17) and (2.21) © 
5 > 7 —aK1 > Anz > 0 and qf > an3 >0 case Eqs.(2.17) and (2.21) (A.21) 


Because (i) in the interval 0 < a < $,sinq@ increases monotonically from 0 to 1 as a increases from 0 to 4; 


and (ii) in the interval 0 < a < 5,0 < sina < 1; and (iii) sin(w — a) = sina for any real number a, one 
concludes that (i) Eq.(A.20) > 


3 
ONGRG ft cin Se SO: (cea) (A.22) 
sin Ap 3 2 
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and (ii) Eq.(A.21) > 


sina,  sin(m — ap) ot al . 
- 1 -)=— : Eq.(A.21 A.2 
sin Oe es > 1 and sin(7) 7) >sinagg >0 (case Eq.( )) (A.23) 
As such, for the case Eq.(A.20), Eqs.(2.16) and (A.22) = (i) 
1 2 
n> > 1.155 (case Eq.(A.20)) (A.24) 
sin AK3 3 


and (ii) 


2 
7 = —= if and only if ag, = age = OK3 = i 
v3 . . (A.25) 
iLe., Qy = A. = 03 = 3 (case Eq.(A.20)) 
On the other hand, for the case Eq.(A.21), Eqs.(2.16) and (A.23) > 
1 2 
—— > V2 = 1.414 > — (case Eq.(A.21)) (A.26) 
sin AK3 V3 


Eqgs.(2.22) and (2.23) now follow directly from Eqs.(A.20)—(A.26). 
To prove Eqs.(2.24) and (2.25) note that both Eqs.(2.20) and (2.21) > 


. > arpa >axn3 >0 for case Eqs.(2.20) or case Eq.(2.21) (A.27) 


By using the fact that, in the interval 0 < a < $,cot a increases monotonically from 0 to +00 as 
a@ decreases from } to 0, (i) Eq.(A.27) > 


cot a3 > cotag2 >0 for case Eqs.(2.20) or case Eq. (2.21) (A.28) 
and (ii) 
lim cota = +00 (A.29) 
aot 


On the other hand, by combining Eqs.(2.19) and (A.28), one has 
2cotan3 > > cotax3 >0 for case Eqs.(2.20) or case Eq.(2.21) (A.30) 


In turn, with the aid of Eqs.(A.27) and (A.29), Eq.(A.30) = both Eqs.(2.24) and (2.25). QED. 


Appendix B 


In this appendix, Eqs.(3.185)—(3.187) will be proved using Eq.(3.139). 
To proceed, note that, with the aid of Eq.(3.50), Eq.(3.139) > 


(a1, A2, 03) = y(a1, 2,7 — a1 — a2) 


def . 9 - 2 - 2 (B.1) 
= g(a1,Q2) = sin* a, + sin“ ag + sin“(a1 + a2), (a1, Qa2,a3) € Da 
where the domain of the function g is 
D(g) 2 {(a1, a2)|a1 > 0,a2 > 0 and a; + a2 < 7} (B.2) 


Given Eqs.(B.1) and (B.2), next we will prove Proposition B-1: 
Proposition B-1: Within its domain D(g), the global maximum of the function g is 2 and it is attained 
at and only at 


(a1, a2) = (1/3, 7/3) (B.3) 
Proof: By using Eq.(B.1) and some trigonometric identities, one has 
0 
= sin(2a1) + sin[2(a1 + a2)] = 2[sin(2a, + a2)](cos a2) (B.4) 
1 


NASA/TM—2017-219534 35 


oo = sin(2a2) + sin[2(a, + a2)) = 2[sin(a; + 2a2)](cos a1) (B.5) 
2 
2 
ua = 4cos(2a, + a2)(cos a2) (B.6) 
Oa 
2 
a = 4cos(a1 + 2a2)(cos a) (B.7) 
005 
and 92 
— 2 cos[2(a1 + a2)] (B.8) 
001002 


By using the last expression in each of Eq.(B.4) and (B.5), it is seen that, within D(g), the extremum 
condition [Ref.12, p. 304] ‘ ‘ 

Dar = Bag 7 (o1102) € D(a) (B.9) 
<> one of the following four possible cases: 
(a) cosa, = cosaz = 0; (b) cosag = 0 and sin(a; + 2a2) = 0; (c) cosa; = 0 and sin(2a; + a2) = 0; and 
(d) sin(2a1 + a2) = sin(a, + 2a2) = 0 

In D(g),0 < a1,a2 < 7. Thus case (a) & ay = ag = 7/2 > a, + QQ = 7. On the other hand, according 
to Eq.(B.2), ay + a2 = 7 => (a1, 02) ¢ D(g). Thus case (a) cannot occur for any (a1,a2) € D(g). 

For any (a1,a2) € D(g),cosaz = 0 & ag = 7/2. On the other hand, with a2 = 4, sin(a, + 2a2) = 
sin(m +a ,) = —sina,. Because sina, > 0 for any (a,,a2) € D(g), case (b) also cannot occur for any 
(a1,a2) € D(g). Similarly one can show that case (c) also cannot occur for any (a1, a2) € D(g). 

Because 2a, + a2 = (a1; + a2) + ay and a, + 2a2 = (a1 + a2) + Qe, it is seen that, for any (a1,a2) € 
D(g),0 < 2a, +a < 2m and 0 < aj + 2a2 < 2m. Thus case (d) © 


2a, tag =a, +20, =7 (B.10) 


for any (a1,a@2) € D(g). Because Eq.(B.10) <= Eq.(B.3), one concludes that, within D(g), the extremum 
condition Eq.(B.9) is met if and only if a1 = ag = §. 
Let a1 = a2 = 3. Then, Eqs.(B.5)-(B.7) > 


0°g ~~ 0Pg 0g 
002 OaR oa Jda,0aq p= oe= ays) ay 
In turn, Eq.(B.11) > 
0?g Og 0" : 
= = = B.12 
daz O03 (G5) Se Aaa am) ( ) 
and 
07g 
Baz ~ —2<0 (a, =a2=7/3) (B.13) 


when a; = ag = 7/3. Because (i) within D(g), the extremum condition Eq.(B.9) is met if and only if 
ay = a2 = 7/3; (ii) Eqs.(B.12) and (B.13) = a relative maximum of g occurs when a1 = a2 = 7/3 [Ref.12, 
p. 312]; (iii) D(g) is an open domain [Ref.13, p.28]; and (iv) 


9 
g(x /3,x/3) = 9 (B.14) 
one concludes that the global maximum of g in D(g) is 9/4 and it is attained at and only at (a,,a2) = 
(x/3, 7/3). QED. 
With the aid of Proposition B-1. Eqs.(B.1) and (B.2) now imply (i) Eq.(3.186), (ii) 


0 < g(a1, a2) < 9/4 if (a1,a2) € D(g) and (a1, a2) 4 (7/3, 7/3) (B.15) 

and (iii) 
lim a1,02) = lim Qa 1,02) = li ; =0 B.16 
ents en 1, 2) ee ee 1,2) OO a ee a2) ( ) 
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Ly 


Q1 
(0,07 “Nz, 0) 
Ly 


Figure 11. Representation of a triangular element in a; — a2 space. 


Moreover, because (i) g is continuous over D(g), and (ii) points (0,0), (0,7) and (7,0) on the a; — ag 
plane (see Fig.11) are limit points [Ref.13, p.28] of D(g) one concludes from Eqs.(B.14)—(B.16) that, given 
any number z with 0 < x < 9/4, there exists an (a2,a2) € D(g) such that g(az,a2) = az. Thus, the range 
of g over D(g) is 0 < g < 9/4. In turn, with the aid of Eqs.(B.1) and (B.2), one concludes that the range of 
y over Dg is defined by Eq.(3.185). 

To prove Eq.(3.187), first note that continuity of g over D(g) coupled with Eqs.(B.14) and (B.15) implies 


that, for any given fixed point (a?,a3) € D(g), (i) g(a?,a$) is a fixed number > 0, and (ii) 


lim | g(a1,a2) = g(a}, a3) > 0, (1,2) € D(g) and (a4, a5) € D(g) (B.17) 


(2,02) (a4 ,02) 


In turn, Eqs.(B.1), (B.2) and (B.17) imply that, as long as (af, a9) is a fixed poin € D(g), it is impossible 
that y > 07 as (a2, a2) > (af, a$) with (a1,a2) € D(g). It becomes possible only if (a?, a) is replaced by 
some special fixed point that lies on the boundary of open domain D(g) so that it lies outside of D(g) and 
yet is a limit point of D(g). 


To search for all the special limit points referred to above, let (i) 


9(Q1, a2) 2 sin? ay + sin? a + sin?(a, +a2) (a1, a2) € D(G) (B.18) 


where, 
D(g) ae {(a1,Q2)|a, > 0,ag >0, and a; + ag < 7} (B.19) 
and (iii) the line segments £1, Lz and L3 depicted in Fig.11 be defined by 


Ly = {(a1,@2)|0 < a, < m and az = 0} (B.20) 
Lo = {(a1,a2)|0 < a2 < m and a; = 0} (B.21) 

and 
Ls © {(a1,a2)| 0 < ay < m and a1 +02 = 7} (B.22) 


Then, (i) one can conclude from Eqs.(B.1), (B.2), (B.18) and (B.19) that the functions of g and @ are identical 
except that the domain D(g) is only a proper subset of D(g); (ii) 


D(g) = D(g) U (D(9))’ (B.23) 


where (D(g))’ denotes the boundary of D(g), i.e., the set of all the limit points of the open domain D(g) 
which lie outside of it, i.e., 
D(g) N(D(g))' = 2 (B.24) 


where, @ denotes the empty set; (iii) By its definition Eq.(B.19), D(g) contains all its own limit points and 
thus it is a closed domain [Ref.13, p.28]; and (iv) with the aid of Eqs.(B.20)—(B.22), the definition of (D(g))’ 
and Fig.11, one concludes that 
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nd 
: LI, N Lz = {(0,0)}, £1 A Ls = {(7,0)}, and L2N Ls = {(0,7)} (B.26) 


Next, let 
s(a@) Wosinta, O<aK<rt (B.27) 
Then, by using Eqs.(B.18)—(B.22), one has 


g(a1, @2) = s(ai), O< ay < wif (a1,Qa2) € Ly (B.28) 

9(a1, @2) = s(az), O< ae < wif (a1, a2) € Ig (B.29) 
and 

g(a1, @2) = s(ai), O< ay < wif (a1, a2) € Ls (B.30) 


On the other hand, Eq.(B.27) and the definitions 


2 
5(q@) a ue and &(a) oe ung O<a<an (B.31) 
da? 
=> 
$(a) = 2sin(2q@) and &(a@) = 4cos(2a), 0<a<a (B.32) 
Eqs.(B.32), in turn, implies that 
Sa) =O0Sa=7/2 if0<a<a7 (B.33) 
and 
§(7/2) =—-4 <0 (B.34) 


By combining Eqs.(B.27), (B.33) and (B.34), one concludes that, (i) no local (relative) minimum of s exists 
in the interval 0 < a < 7, (ii) the global maximal value of s in the interval 0 < a < 7 is s(m/2) = 2, and it 
is attained if and only if a = 7/2, ie., 


s(77/2) =2>s(a) >0if0<a<7/2or7/2<a<7 (B.35) 


and (iii) 

s(0) = s(7) =0< s(a), 0O<a<m7 (B.36) 
With the aid of the results presented in the above items (i)—(iii), Eqs.(B.25) and (B.28)—Eqs.(B.30) imply 
that, over the boundary (D(g))’ of the open domain D(g), (i) the maximal value of g is 2, 
and it is attained at and only at (a1,Q@2) = (7/2,0), or (@1,Q@2) = (0,7/2), or (a1,Qa2) = 
(2 /2,7/2); and (ii) the minimal value of g is 0, and it is attained at and only at (a1,a2) = 
(0,0) or (@1,a2) = (7,0) or (a1,a2) = (0,7). Because (a) the function g is continuous over its 
domain D(g), and (b) the points (0,0), (7,0) and (0,7) on the a; — a2 plane are the only limit points of the 
open domain D(g) which lie in (D(g))’ and yet the function 9 attains its minimal value 0 at these points, 
Eq.(3.187) now follows directly from Eqs.(B.1), (B.2), (B.15)-(B.19) and (B.23). Note that the fact that the 
maximal value of g attained in (D(g))’ is 2, is consistent with Eq.(3.186). 
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